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Abstract—Postflexion and transfor- 
mation stages of the Diaphus ostenfeldi 
in the southwestern Pacific Ocean are 
described for the first time. Identifica- 
tion to the species level was based on 
tracing characters from fully formed, 
large metamorphic specimens back- 
ward to smaller specimens. It was 
found that larval D. ostenfeldi =9.5 mm 
standard body length (SL) can be char- 
acterized by a unique character within 
the genus Diaphus, namely the pres- 
ence of melanophores at the dorsal fin 
base and between branchiostegal rays. 
Metamorphosis in larvae of this species 
occurs at a much longer body length 
(~20 mm SL) than in larvae of other 
species of this genus. 
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Diaphus is the most speciose genus of 
the mesopelagic family Myctophidae 
and plays a key role in oceanic food 
webs. Larvae of Diaphus are poorly 
known. Of the 80 species in this genus 
(Fricke et al., 2020), the early stages of 
development are known for only 10 spe- 
cies. These species are the D. agas- 
sizit (Pertseva-Ostroumova, 1975), 
Garman’s lanterfish (D. garmani) (Sassa 
et al., 2003), D. hudsoni (Olivar, 1987), 
short-headed lantern fish (D. brachy- 
cephalus or D. richardsoni) (Olivar 
and Beckley, 1995), D. diadematus 
(Olivar and Beckley, 1995), D. mollis 
(Olivar and Beckley, 1995), D. meto- 
poclampus (Sparta, 1952), D. pacificus 
(Moser and Ahlstrom, 1974), Califor- 
nia headlightfish (D. theta) (Moser and 
Ahlstrom, 1996), and Malayan lantern- 
fish (D. malayanus) (Tsokur, 1975). 
The early stages of development are 
known in only 2 species of the 9 species 
of Diaphus that have been observed 
in the southwestern Pacific Ocean: 
D. hudsoni and D. mollis (McGinnis, 
1982; Bekker and Evseenko, 1986). 
Diaphus ostenfeldi is circumglobal, 
is distributed in the Subtropical Con- 
vergence, a frontal zone between sub- 
antarctic and tropical water masses, 
and is considered to be a convergence 
(Krefft, 1974), notal or subantarctic 
(Parin et al., 1974), or transitional 


(McGinnis, 1982; Bekker and Evseenko, 
1986) species. This species has been 
observed in the central eastern Pacific 
Ocean (McGinnis, 1982) and in the 
southwestern Pacific Ocean (Bekker 
and Evseenko, 1986), and its range is 
known to extend to about 23-24°S in 
the region of the Benguela Current 
(Trunov, 1968; Hulley, 1972). It has been 
found at depths as great as about 120 m 
below the sea surface (Bekker, 1983). 
Juveniles have been taken at depths 
less than 100 m; larger specimens have 
been found at depths 160 m below the 
sea surface (Hulley, 1981; McGinnis, 
1982). Little is known of the ecology 
of this species. Specimens caught in 
the South Atlantic Ocean were limited 
between the 10°C and 15°C isotherms 
at a depth of 200 m; however, it has 
been assumed that the species can live 
at shallower depths, in suitable tem- 
peratures (Krefft, 1974; Hulley, 1981). 

Herein, we present larval descriptions 
for D. ostenfeldi from the southwestern 
Pacific Ocean based on a developmen- 
tal series of larvae from the postflexion 
stage to the transformation stage. 


Materials and methods 


Larvae of D. ostenfeldi were caught at 
4 stations (stations [st.] 3049, 3052, 
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3054, and 3069) during the 34th expedition of the RV Dmi- 
try Mendeleev to the southwestern Pacific Ocean in Jan- 
uary—February 1985 (Fig. 1). All specimens were caught 
by using an Isaacs-Kidd mid-water trawl with the Samy- 
shev-Aseev modification in the depth layers of 0-100 m, 
0-200 m, and 0-600 m in the area between 35°S and 57°S. 
This trawl is 25 m long and has a 5-mm mesh net without 
nodes and with a terminal insertion of 500-ym nylon mesh 
and a 6-m” mouth area. The specimens were preserved and 
stored in 10% formaldehyde. We captured 62 specimens 
that were 8.6—20.5 mm standard body length (SL); larvae of 
smaller sizes were absent in samples of trawl hauls. 
Information for the following features are included in 
the descriptions of larval stages in this paper (Moser, 
1996): SL, head length (HL), body depth at pectoral fin 
base (BD), preanal length (PAL), predorsal length (PDL), 
snout length (SnL), horizontal diameter of the eye (ED), 
number of rays in dorsal fin (D), number of rays in anal 
fin (A), number of rays in pectoral fin (P,), number of rays 
in pelvic fin (P,), number of rays in caudal fin (C), number 
of gill rakers (GR) on the upper and lower parts of the 
first gill arch, and number of vertebrae (V). Descriptions 
include patterns of the following photophores (Bekker, 
1983): branchiostegal (Br), opercular (Op), pectolateral 
(PLO), pectoventral (PVO), pectoral (PO), ventrolateral 
(VLO), ventral (VO), supraanal (SAO), anterior anal 


(AOa), posterior anal (AOp), and posterolateral (Pol). Also 
included in descriptions are the following orbital glands: 
dorsonasal (Dn), anteorbital (Ant), suborbital (So), and 
ventronasal (Vn). 


Results 
Identification 


The identification of large specimens (>20 mm SL) as 
D. ostenfeldi was based on the following features: the num- 
ber of gill rakers on the first gill arch of 25-26, So absent, 
Ant present, Dn and Vn connected but not extended ante- 
riorly from the olfactory rosette, no enlarged teeth in the 
lower jaw inner row, head depth approximately equal to its 
length, VLO about midway between lateral line and base of 
the ventral fin. The identification of specimens smaller 
than 20 mm SL was based on tracing characters from fully 
formed, large metamorphic specimens backward to smaller 
specimens. 


Body proportions 


The larvae have a moderately deep body (BD 24—28% SL). 
The head is relatively small (HL 25-30% SL) but deep with 
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Figure 1 


Schematic map of the stations (open circles) where larval (open squares) and adult (black 
circles) Diaphus ostenfeldi were sampled during the 34th cruise of the RV Dmitry Men- 
deleev in the southwestern Pacific Ocean between January and February 1985. Dashed 
lines indicate the positions of biogeographic borders (Becker and Evseenko, 1986) for the 
central water mass (C), periphery of the central water mass (P), subtropical frontal zone 
(STFZ), notal zone (N), Antarctic Polar Frontal Zone (APFZ), and Antarctic (A). Station 
numbers are given next to station markers. 
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a steep profile; the depth of the head at the corner of the 
mouth is almost equal to HL. The jaws are long and con- 
tinue beyond the vertical of the eye posterior edge; both 
jaws have small teeth. Jaw length increases with larval 
size. The snout is short (SnL 19-22% HL), the eye is slightly 
oval at sizes up to about 16 mm SL (ED 25-30% HL). The 
anus opens slightly behind the midbody (PAL 57-61% SL). 
The dorsal fin begins almost at the middle of the body (PDL 
4448% SL). The posterior margin of the pectoral fin rays 
reaches the vertical of the ventral fin base, and the tips of 
the ventral fin rays reach the anal fin origin. The propor- 
tions are generally retained with increases in larval size. 


CMOOM oe 


Meristic characters 


D: 16-17, A: 16, P,: 11, P,: 8, C: 6—-10+9-6(5), GR: 
(8)9+14+15—-17=25-26, V: 37-38. Fin rays are formed in all 
fins at 8.6 mm SL, and larvae have a definitive number of 
gill rakers at 12.0 mm SL. Larvae have 37-38 myomeres. 


Pigmentation 
The 8.6-mm-SL larva (Fig. 2A) has one melanophore dor- 


sally on the peritoneum near the terminus of the gut, one 
melanophore behind the anal fin base (these melanophores 


Ny WN WSS — —. 


Figure 2 


Illustrations of 5 larval and transformation specimens of Diaphus ostenfeldi caught in the southwestern 
Pacific Ocean in 1985: (A) 8.6 mm standard body length (SL), (B) 9.5 mm SL, (C) 10.7 mm SL, (D) 10.7 mm 
SL (dorsal view), (E) 10.7 mm SL (ventral view), (F) 15.8 mm SL, and (G) 20.5 mm SL. 
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are present in all specimens and noticeable to a size of 
20.5 mm SL), and one melanophore above the posterior 
part of the anal fin base (approximately at the level of the 
12th pterygiophore). Several melanophores are present 
proximally on the lower caudal fin rays. Internal melano- 
phores are visible above the midbrain and above the otic 
capsules; melanophores are scattered on the dorsolateral 
part of the peritoneum. In the 9.5-mm-SL larva (Fig. 2B), 
melanophores are added at the dorsal fin base and at the 
base of some upper caudal fin rays. An unpaired internal 
melanophore is visible on the anteroventral part of the 
liver on the level of the pectoral fin base; melanophores 
appear above the hindbrain, on the opercle, and on the 
membrane between the branchiostegal rays. Pigmenta- 
tion becomes more intense by a size of 10.7 mm SL (Fig. 
2, C—E) and is added ventrally at some myosepta on each 
side along the anal fin base. Melanophores at the dor- 
sal fin base and above the anal fin base are not located 
symmetrically on each side. This pigmentation pattern 
persists up to a size of 20.5 mm SL (Fig. 2, F and G). Ven- 
trally, almost all individuals larger than 10.7 mm SL have 
a row of 2—4 melanophores: approximately at the level of 
the future PO,, PO;, VO,, and VO;. Pigmentation of the 
caudal fin persists to a size of 20.0 mm SL. 


Photophore formation 


The 8.6-mm-SL larva has Br,, PO,, PO;, and VO, begin- 
ning to form. At a size of 11.0 mm SL, Br, and PO, ; are 
added. Larvae at a size of 13.0 mm SL have Br, 3, Op, », 
PO,_;, PVO,_5, VO,_;, and AOa, 7. At a size of 16.0 mm 
SL, formation of PLO, VLO, SAO, _3, Pol, and AOp,_; ends 
(Fig. 2). By a size of 20.0 mm SL, specimens have a defini- 
tive photophore pattern, but formation of the orbital lumi- 
nous glands is not yet complete: So absent, the dark tissue 
at the future Ant is triangular, roundish Dn connected 
with the dark tissue of the future Vn and Ant, and Vn 
reaches the anteroventral part of the eye (Fig. 3). 


Discussion 


The larvae described herein belong to a moderately deep 
body morphotype (Moser and Ahlstrom, 1972, 1974; Moser 
et al., 1984). They have a few melanophores on the pos- 
tanal ventral midline (Fig. 2) and have no So photophore 
(Fig. 3). In total, 9 species of Diaphus can be found in our 
study area in the southwestern Pacific Ocean: D. ander- 
seni, D. danae, D. effulgens, D. hudsoni, small lantern fish 
(D. meadi), D. mollis, Taaning’s lantern fish (D. termophi- 
lus), D. parri, and D. ostenfeldi (McGinnis, 1982; Bekker 
and Evseenko, 1986). Larvae are known for 2 of these spe- 
cies, D. hudsoni and D. mollis (Olivar, 1987; Olivar and 
Beckley, 1995), and they belong to a different morphotype 
that is moderately slender. In addition, they clearly differ 
from our larvae by the absence of pigmentation at the dor- 
sal fin base and at the branchiostegal rays. Among the 
remaining 7 species, only 3 species, D. effulgens, D. danae, 
and D. ostenfeldi, have no So photophore as we observed in 


Figure 3 


Illustration of the locations of the orbital luminous organs 
in a transformation specimen of Diaphus ostenfeldi that 
has a standard body length of 20.5 mm and was caught 
in the southwestern Pacific Ocean in 1985. Locations of 
the following organs are shown: dorsonasal (Dn), ante- 
orbital (Ant), ventronasal (Vn), branchiostegal (Br), and 
opercular (Op). 


our transformation specimens. One of these species, 
D. effulgens, has fewer gill rakers on the first gill arch than 
we observed in our specimens (17-22 versus 24-25); 
another of those species, D. danae, has more gill rakers on 
the first gill arch (27-30) (Nafpaktitis, 1978; Bekker, 
1983). For these reasons, the Diaphus larvae in the devel- 
opmental series presented herein were identified as lar- 
vae of D. ostenfeldi. 


Comparative notes 


There are several common features in the pigmentation 
between our Diaphus larvae and those described in earlier 
reports (Sparta, 1952; Moser and Ahlstrom, 1974, 1996; 
Pertseva-Ostroumova, 1975; Tsokur, 1975; Olivar, 1987; 
Olivar and Beckley, 1995; Sassa et al., 2003): melano- 
phores behind the cleithrum symphysis, at the peritoneum 
near the terminus of the gut, behind the anal fin base, and 
proximally on the lower caudal fin rays. At the same time, 
larval D. ostenfeldi can be distinguished from all known 
Diaphus larvae by the presence of pigment on the dorsal 
fin base and the branchiostegal rays at lengths 29.5 mm 
SL. One more difference between the larvae of D. ostenfeldi 
and the other Diaphus species is that parietal and perito- 
neal pigmentation appears earlier in development in the 
D. ostenfeldi than in its congeners. In addition, the trans- 
formation of larvae described in this article was still not 
complete even at a length of 20.5 mm SL; our larvae have 
a full complex of photophores typical for adults, but the 
formation of orbital glands was not yet completed. In all 
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Diaphus larvae described earlier by others (Sparta, 1952; 
Moser and Ahlstrom, 1974, 1996; Pertseva-Ostroumova, 
1975; Tsokur, 1975; Olivar, 1987; Olivar and Beckley, 1995; 
Sassa et al., 2003), photophores are formed by a size of 
about 11-15 mm SL. This length may correlate with the 
size of adults. The Diaphus ostenfeldi attains a length of 
about 120 mm SL (Hulley, 1981), but the sizes of other spe- 
cies whose early stages are known do not exceed 90 mm 
SL (Hulley, 1981; Bekker, 1983). 


Distribution 


In the southwestern Pacific Ocean between 35°S and 57°S, 
adults of this species were encountered at the 18 stations 
sampled during the 34th expedition of the RV Dmitry 
Mendeleev in different water masses, from the central 
masses in the north of the study area across the subtrop- 
ical frontal zone and the notal zone of the South Pacific 
Ocean to the Antarctic Polar Frontal Zone in the south 
of the study area (Bekker and Evseenko, 1986). Larvae of 
D. ostenfeldi were found at 4 stations in a narrow range of 
latitudes between 35°S and 39°S in the subtropical frontal 
zone. This area was located on the western periphery of 
the vast anticyclonic South Pacific Gyre. 
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Abstract—The stone scorpionfish (Scor- 
paena mystes) is common throughout 
the Gulf of California in Mexico. Because 
little information is available regard- 
ing the basic biological characteristics 
of this species and because commercial 
demand for it is increasing, the need 
for studies of population dynamics is 
urgent. We sampled 233 stone scorpi- 
onfish (117 males and 116 females) in 
the Gulf of California from May 2015 
through April 2016 to estimate age and 
individual growth. Total lengths (TL) 
ranged between 15.3 and 35.5 cm for 
males and between 19.0 and 44.5 cm 
for females. Results from analysis of 
the edge types of otoliths (opaque and 
translucent) indicate an annual period- 
icity in the formation of annuli in sagit- 
tae. Ages of stone scorpionfish ranged 
from 2 to 10 years for males and from 
3 to 15 years for females. Estimates 
for the parameters of the von Berta- 
lanffy growth function were 34.76 and 
48.22 cm TL for asymptotic length, 
0.202/year and 0.135/year for growth 
coefficient, and -1.28 and —0.80 years 
for theoretical age at length zero for 
males and females, respectively. This 
study resulted in the first estimates 
of age and individual growth of stone 
scorpionfish, information that is essen- 
tial for evaluation and management of 
this important fishery resource. 
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The stone scorpionfish (Scorpaena mys- 
tes) is commonly found throughout the 
Gulf of California in Mexico, with a dis- 
tribution in the eastern Pacific Ocean 
that extends from southern California to 
Chile, including offshore islands (Butler 
et al., 2012). This demersal species occurs 
in shallow waters and seaweed-covered 
reefs as well as in open, sandy areas to 
depths of 85 m. It is considered an oppor- 
tunistic predator and feeds mainly on 
crustaceans and small fish (Butler et al., 
2012). Adults reach sizes up to 51 cm in 
total length (TL) (Robertson and Allen, 
2015). 

The meat of stone scorpionfish is 
considered to be of excellent quality 
and has high commercial value; as 
a result, this species has become an 
important fishery resource (Jakes-Cota 
et al., 2017). Historically, however, 
this species has not been sought after 
commercially because venom glands 
located along the spines of the fins 
can cause painful injuries and severe 
respiratory disorders (Poss, 1995); it 
has long been a secondary priority spe- 
cies in the artisanal fishery in the Gulf 
of California. Since 2012, decreases in 
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the catches of other species of high 
commercial value, such as the Pacific 
red snapper (Lutjanus peru) and yel- 
lowtail jack (Seriola lalandi), caused 
an increase in price and capture vol- 
ume of the stone scorpionfish, driving 
this species into the commercial sec- 
tor as a complementary, high-quality 
product. 

Given the little amount of data avail- 
able on the basic biological characteris- 
tics of this species and its recent 
commercial exploitation in the region 
of the Gulf of California, there is an 
urgent need to document its basic life 
history (e.g., age, growth, mortality, 
and reproduction) and its population 
dynamics and to create a management 
protocol for this fishery resource. Age 
and individual growth should be 
included as key elements in every 
investigation focused on the rational 
exploitation of fishing resources 
(Beamish and Fournier, 1981). Prior to 
this study, the age and individual 
growth characteristics of stone scorpi- 
onfish were not available; therefore, 
the aims of our study were 1) to esti- 
mate age by using otoliths, 2) to 
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estimate individual growth parameters, 
and 3) to test for significant differences 
in growth trajectories between males 
and females. 


Materials and methods 


Stone scorpionfish were collected once 
a month (21—22 individuals per month) 
in the artisanal fishery from May 2015 
through April 2016. The fish were 
taken by spearing while using semiau- 
tonomous diving equipment (hookah), 
at depths up to 30 m, at different fish- 
ing sites along 20 km of coastline south 
of the harbor of Santa Rosalia in Baja 
California Sur, Mexico. The standard 
length (SL) and TL were measured to 
the nearest millimeter, and the weight 
(W) was recorded to the nearest 0.1 g 
for each individual. Sex was deter- 
mined by using macroscopic exam- 
ination of the gonads, and the exact 
binomial test was used to test differ- 
ences in sex ratio. Pairs of sagittae were 
extracted, cleaned, and stored dry in 
plastic vials in the laboratory, and only the right otoliths 
were used for the analyses described in the rest of this 
section. 

To describe the length—weight relationship, we used this 
equation: 


W =aTL’, (1) 


where a = the intercept of the regression line; and 
b = the allometric coefficient, the slope of the regres- 
sion line. 


To describe the relationship independently for males and 
females with log-transformed length and weight data, we 
used this equation: 


log W =a + blog TL. (2) 


The a and 6 parameters were estimated by using a sim- 
ple linear regression analysis. A Student’s t-test was 
performed to evaluate whether the b-value was signifi- 
cantly different from 3, a value that indicates an isometric 
growth pattern (i.e., all body parts grow at approximately 
the same rate, and a fish has an unchanging body form 
throughout development). An analysis of covariance was 
used to test differences in b-values between males and 
females. 

The otolith radius (OR) (Fig. 1) was measured under 
a stereoscope (Stemi SV11', Carl Zeiss AG, Oberkochen, 
Germany) at 3.2x magnification by using the imaging 
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Figure 1 


Photograph showing the annular pattern in the sagittae of a 4-year-old male 
stone scorpionfish (Scorpaena mystes) caught in May 2015 in the Gulf of 
California in Mexico. Black dots indicate annuli. The black arrow indicates the 
measurement of the otolith radius. The scale bar is equal to 1 mm. 


software AxioVision, vers. 4.6 (Carl Zeiss AG). The ste- 
reoscope was attached to a digital camera (AxioCam 
MRc 5, Carl Zeiss AG). The relationship between TL and 
OR was described by using a simple linear regression 
analysis. 

Right otoliths were immersed in ethanol with the con- 
cave side down, and 2 independent readers counted the 
annuli once with a stereoscope, using reflected light 
against a dark background without prior knowledge of 
fish length or weight. Under reflected light, the nucleus 
(the core of the otolith) and the opaque bands appeared 
as light rings and the translucent bands appeared as dark 
rings (e.g., La Mesa et al., 2005, 2010; Bilgin and Celik, 
2009). The combination of an opaque band followed by 
a translucent band was considered to be an annulus, as 
has been done for other scorpionfish species (e.g., Massuti 
et al., 2000; La Mesa et al., 2005, 2010; Bilgin and Celik, 
2009). Annuli were counted from the nucleus toward the 
tip of the rostrum (the anterior projection of the otolith) 
along the same line (Fig. 1). The precision between readers 
was evaluated by using the mean coefficient of variation 
(CV) and average percent error (APE). Values of 5.5% for 
APE and of 7.6% for CV have been considered adequate in 
many aging studies (Campana, 2001). 

To validate the periodicity of the formation of opaque 
and translucent bands and to assign a unit of time (e.g., 
years) to the counted annuli, edge type analysis was 
used. For this analysis, the edge type of each otolith was 
recorded as either opaque or translucent, and, for each 
month, the percentages of opaque edges and of translu- 
cent edges in the total number of otoliths analyzed was 
plotted. 
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The von Bertalanffy growth function was fitted to 
length-at-age data to estimate individual growth: 


L, = L(t = ee) (3) 


where L, = the length (TL, in centimeters) at age ¢; 
L,, = the asymptotic length; 
k = the growth coefficient; and 
ty = the theoretical age at length zero. 


Additionally, a likelihood ratio test (Kimura, 1980), also 
known as the likelihood ratio chi-square test, was applied 
to test for differences between males and females in 
parameters of individual growth. 


Results 


During the study period, 233 stone scorpionfish (117 males 
and 116 females) were sampled. No significant differ- 
ence was found in the ratio of males to females when 
applying a 1:1 exact binomial test (P>0.05). However, the 
number of males was lower than the number of females 
for the large fish (>30 cm TL). Overall, 70% of the 
sampled individuals ranged in size between 20 and 
30 cm TL. The lengths of the males ranged from 15.3 to 
35.5 cm TL (12.5—-28.9 cm SL), with a mean length of 
25.9 cm TL (SE 3.60). Comparatively, the lengths of the 
females ranged from 19.0 to 44.5 cm TL (15.5-36.5 cm 
SL), with a mean length of 27.8 cm TL (SE 5.17). Weights 
ranged from 64 to 1033 g for males, with a mean weight 
of 406 g (SE 181), and from 119 to 1870 g for females, 
with a mean weight of 545 g (SE 385). The analysis of 
variance revealed significant differences in mean length 
(Fy.931)=10.67, P<0.05) and mean weight (F(, 93;)=12.57, 
P<0.05) between males and females. Females were longer 
and heavier than males. 

The length—weight relationship (Table 1) was not signif- 
icantly different between males and females (F(; 999)=1.27, 
P>0.05). The b-value of the length—weight relationship of 
males and females was not significantly different from 3 
(P>0.05), indicating an isometric growth pattern. 


Otolith radius ranged between 2.5 and 5.6 mm. The lin- 
ear regression model developed to predict the relationship 
between OR and TL was as follows: 


TL = -3.53 + 3.910R. (4) 


The linear relationship between both variables was signif- 
icant (coefficient of determination [r7]=0.79, P<0.05). 

The highest percentage of otoliths with opaque edges 
was recorded in August (81%), followed by the percentage 
in December (65%), and the highest percentages of otoliths 
with translucent edges were observed in February (80%) 
and June (68%). These results indicate that the periodicity 
of annulus formation is annual and that the formation of 
a new annulus begins with an opaque band in the months 
of August-September and ends with a translucent band in 
the months of May—July. 

The detections of annuli were confirmed and annuli 
could be counted in 98% of the 233 otoliths analyzed. 
For the remaining 2% of otoliths, the annuli could not be 
counted because they were not visible. The ages of individ- 
ual specimens ranged between 2 and 10 years for males 
and between 3 and 15 years for females. The most fre- 
quent ages ranged from 5 to 7 years (63%) for males and 
from 4 to 6 years (74%) for females. The values of CV and 
APE were 6.50 and 4.59, respectively. 

The parameters of the von Bertalanffy growth function 
for males and females and the growth curves that were 
fitted to length-at-age data are shown in Figure 2. Results 
of the likelihood ratio test indicate significant differences 
in parameters of the von Bertalanffy growth function 
between males and females (y7=1.89, P<0.05); females 
grew slower (indicated by a lower k value) yet reached a 
greater maximum average length (indicated by a greater L,,.) 
than males. 


Discussion 


In this study, the maximum length recorded for stone scor- 
pionfish was 44.5 cm TL, which is comparable to the maxi- 
mum length of 51 cm TL previously reported for this species 
(Robertson and Allen, 2015). In addition, individuals larger 


Table 1 


Estimates of the parameters from regression analysis of the length—weight relationship for 
male and female stone scorpionfish (Scorpaena mystes) caught in the Gulf of California in 
Mexico during May 2015—April 2016. The parameters are the intercept (a) and slope (6) of the 
regression line. r’=coefficient of determination; ClI=confidence interval; SE=standard error of 


the mean. 


Sex Mean SE 


Males 0.0205 
Females 0.0159 


0.29136 
0.18687 


95% CI 


0.0115-0.0365 
0.0109-0.0229 


b 


95% CI r? 


Mean SE 


3.01 0.08975 
3.09 0.05696 


2.84-3.19 0.90 
2.98-3.21 0.93 
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—® Males: L, = 34.76(1 — e-0-202(-(-1.28))) 


—@— Females: L; = 48.22(1 —e~? 195-08) 


Results of the edge type analysis indi- 
cate that the periodicity of annulus for- 
mation in stone scorpionfish has an 
annual pattern. The highest percentages 
of opaque edges among all otoliths ana- 
lyzed coincided with the highest condi- 
tion factor values reported by Jakes-Cota 
et al. (2017) for stone scorpionfish for 
the months of August and September. 
This finding indicates that conditions for 
growth (e.g., high temperatures and high 
food availability) are favorable in our 
study area (Le Cren, 1951). In contrast, 
the highest percentages of translucent 
edges among all otoliths analyzed coin- 
cided with the lowest condition factor val- 
ues reported by Jakes-Cota et al. (2017) 
for the months of February—June. This 
finding indicates conditions (e.g., low 
temperatures and scarce food) are unfa- 
vorable in our study area and may limit 
growth (Le Cren, 1951). The annual peri- 


T 
10 
Age (years) 


Figure 2 


Growth curves from the von Bertalanffy growth function fitted to length- 
at-age data for male (black dots) and female (gray dots) stone scorpionfish 
(Scorpaena mystes) caught in the Gulf of California in Mexico during May 


2015—April 2016. L,=total length at age t. 


than 30 cm TL were scarce in our sampling, possibly a 
consequence of cumulative mortality in which larger and 
older individuals are removed by fishing or die naturally 
(Csirke, 1980). 

The characteristics of the length—weight relationship 
indicate an isometric growth pattern for both males and 
females. Similar results were reported for this species 
even when males and females were not analyzed sepa- 
rately (Jakes-Cota et al., 2017). In contrast to these find- 
ings, it has been reported that other species of the genus 
Scorpaena, such as the black scorpionfish (S. porcus) (Bilgin 
and Celik, 2009; Demirhan and Can, 2009), Madeira rock- 
fish (S. maderensis) (La Mesa et al., 2005), and small red 
scorpionfish (S. notata) (Neves et al., 2015), commonly 
have positive allometric growth. However, as in our study 
of the stone scorpionfish and in studies of other scorpion- 
fish species, such as the blackbelly rosefish (Helicolenus 
dactylopterus) (Massuti et al., 2000), results also indicate 
isometric growth. 

The relationship between OR and TL of stone scorpion- 
fish in our study was linear and statistically significant, 
indicating that these measurements have a relation of 
proportionality and that the OR may be a predictor of fish 
length. This relationship has been reported for other species 
of Scorpaena, such as the Madeira rockfish (La Mesa et al., 
2005) and small red scorpionfish (Scarcella et al., 2011). 


odicity of annulus formation in otoliths 
has been reported for other scorpionfish 
species, such as the black scorpionfish 
(Bilgin and Celik, 2009; Demirhan and 
Can, 2009; La Mesa et al., 2010), small 
red scorpionfish (Scarcella et al., 2011; 
Neves et al., 2015), Madeira rockfish 
(La Mesa et al., 2005), and blackbelly 
rosefish (Massuti et al., 2000). 

Annuli were apparent in the otoliths 
of stone scorpionfish used in this study. 
In general, translucent bands are slightly wider than the 
opaque bands, and as the number of annuli increased, 
the widths of both bands decreased. Both bands were 
clearly observed and distinguished, and it was relatively 
simple to count them, as evidenced by the low values of 
CV and APE; these low values indicate a high concor- 
dance between counts performed independently by read- 
ers and reflect consistency between observations. In this 
study, we found a significant relationship between the 
OR and length of fish, an annual pattern of annulus for- 
mation in otoliths, and a high precision in the counting of 
annuli between technicians. These results indicate that 
the otolith is a structure that can be used reliably in age 
estimation of stone scorpionfish. 

This study is the first to report age and individual 
growth characteristics for stone scorpionfish. We found 
males up to 10 years old and females up to 15 years old, 
and we determined that individuals 4—7 years old com- 
posed the most abundant age class in our sample. 

Individual growth characteristics were significantly 
different between males and females. In addition, our 
results indicate that females are larger, heavier, and 
older than males. This pattern of distinctive growth 
between males and females has been reported previ- 
ously for the black scorpionfish (Bilgin and Celik, 2009). 
However, it is also common that males reach lengths, 
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ages, and L., that are greater than those of females, as 
has been observed for other species, such as the Madeira 
rockfish (La Mesa et al., 2005), small red scorpionfish 
(Neves et al., 2015), and blackbelly rosefish (Massuti 
et al., 2000). Results from our study indicate that males 
and females take approximately 10 years to reach a size 
of 30 and 35 cm TL, respectively, indicating that the 
stone scorpionfish is a slow-growing species, like most 
scorpionfish species. 

The estimates of ages and individual growth rates 
reported from this study improve our knowledge of the 
stone scorpionfish. Information from this study can be 
used for stock assessment and management of a species of 
emerging fishing potential in the artisanal fishery of the 
Gulf of California. 
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Abstract—Results from published age 
and growth models for bonnetheads 
(Sphyrna tiburo) indicate significant 
differences in life history between pop- 
ulations in the eastern Gulf of Mexico 
(GOM) and those in estuarine waters of 
the Atlantic coast of the southeastern 
United States (hereafter referred to as 
the Atlantic region). An age-indepen- 
dent model, GROTAG, was used with 
region-specific tag-recapture data to 
generate estimates of von Bertalanffy 
growth parameters and growth rates 
for sharks in each of these regions, and 
these estimates were compared with 
age-based life history characteristics. 
Results from the GROTAG model indi- 
cate that female bonnetheads in the 
GOM initially grew faster and attained 
a smaller maximum size than females in 
the Atlantic region. The final GROTAG 
model for females in the Atlantic region 
produced estimates of von Bertalanffy 
parameters and growth rates similar to 
those produced by the age-based growth 
model. For the population in the GOM, 
GROTAG model results indicate that 
growth rates were slower and aver- 
age maximum size and longevity were 
greater than those from age-based mod- 
els. Although models for males were 
generated with tag-recapture data, 
large 95% confidence intervals hin- 
dered comparisons. For both sexes and 
regions, calculated maximum longevity 
and age at 50% maturity are larger than 
published estimates, indicating that age 
underestimation may have occurred in 
both age and growth studies, with sig- 
nificant differences in estimates of life 
history characteristics for bonnetheads 
in the GOM. 
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Growth of wild fish is primarily esti- 
mated by using 3 types of models: 
age-based models in which length and 
estimated age are known; length-based 
models that use modal progressions 
in length-frequency data; and models 
that use length at capture, known time 
at liberty, and length at recapture from 
tag-recapture data to model growth 
trajectories. Length-based models are 
useful; however, they are suitable only 
for estimating growth in younger age 
classes because substantial overlap 
in lengths of fish can occur in older 
cohorts (Campana, 2001; Natanson 
et al., 2018a). Methods based on the 
use of age and tag-recapture data 
often characterize growth of fish by 
using the von Bertalanffy growth func- 
tion (VBGF) (Ricker, 1975; Francis, 
1988a). Growth of fish is most com- 
monly estimated by using age-based 
models; however, models developed 
with tag-recapture data are often used 
as alternative models for large, long- 
lived species that may not be able to 
tolerate high rates of fishing mortality 


4 Northeast Fisheries Science Center 
National Marine Fisheries Service, NOAA 
28 Tarzwell Drive 
Narragansett, Rhode Island 02882 


> Mississippi Laboratories 
Southeast Fisheries Science Center 
National Marine Fisheries Service, NOAA 
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or for species for which age estimation 
is difficult or imprecise (Cailliet et al., 
1992; Treble et al., 2008). 

All methods of modeling growth 
have limitations and potential biases. 
Age estimation can be imprecise and 
is difficult to validate for many species 
(Beamish and McFarlane, 1983), espe- 
cially elasmobranchs (Goldman, 2004) 
and deepwater species (Rigby et al., 
2014). Validation of age estimates is 
recommended for all age and growth 
studies; however, validation of all age 
classes is difficult and rarely achieved 
(Cailliet et al., 2006). Although val- 
idation of ages has been completed 
for several species of elasmobranchs 
(Campana et al., 2002; Ardizzone et al., 
2006; Kneebone et al., 2008), there is 
a growing body of evidence that indi- 
cates that age underestimation com- 
monly occurs for long-lived species, 
such as sharks (e.g., Francis et al., 
2007; Passerotti et al., 2014; Harry, 
2018; Natanson et al., 2018b), as well 
as recent evidence that age underesti- 
mation may also occur in sharks with 
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intermediate longevities (Frazier et al., 2014). For these 
reasons, modeling of growth with tag-recapture data and 
newer methods of modeling that combine tag-recapture 
data with age estimates have been recommended (Eveson 
et al., 2007; Aires-da-Silva et al., 2015; Francis et al., 2016; 
Natanson and Deacy, 2019). 

Growth models that use tag-recapture data also have 
limitations and biases. To accurately represent growth 
in a population, recapture data should include individ- 
uals throughout the range of lengths found in the pop- 
ulation. However, this goal is rarely achieved with most 
data sets. Further, growth of fish may be affected by tag- 
ging (Kelly and Barker, 1963; Saunders and Allen, 1967), 
including elasmobranch growth (Gruber, 1982; Parsons, 
1987; Davenport and Stevens, 1988; Kalish and John- 
ston, 2001). Measurement error in body length of tagged 
and recaptured fish can also bias models; however, some 
methods of modeling data can incorporate measurement 
error as a parameter (Francis, 1988b). When sufficient 
tag-recapture data are available, models can be used to 
verify or compare growth rates and estimates of longev- 
ity (maximum age) and age at maturity with those from 
traditional age-based models (e.g., Natanson et al., 2002; 
Welsford and Lyle, 2005; Natanson and Deacy, 2019). One 
advantage of models based on tag-recapture data is that, 
if recapture data spans multiple seasons, information 
about growth variability and seasonal changes in growth 
can be determined (Francis, 1988b). Such information 
can be especially useful when comparing region-specific 
growth rates between or within populations. 

The bonnethead (Sphyrna tiburo) is a relatively small 
shark species, with individuals reaching a maximum 
size of 150 cm total length, that is commonly found in 
the coastal and estuarine waters of the western North 
Atlantic Ocean from North Carolina to southern Brazil, 
including the Gulf of Mexico (GOM) and the Caribbean 
Sea (Compagno, 1984). Significant differences in life 
history characteristics exist between bonnetheads cap- 
tured off the Atlantic coast of the southeastern United 
States (hereafter referred to as the Atlantic region) 
and those caught in the eastern GOM (Frazier et al., 
2014). Regional variation in life history between popu- 
lations in the Atlantic region and in the GOM has been 
found for other coastal shark species, including the 
Atlantic sharpnose shark (Rhizoprionodon terraenovae) 
(Carlson and Loefer!), blacknose shark (Carcharhinus 
acronotus) (Driggers et al., 2004), and finetooth shark 
(C. isodon) (Drymon et al., 2006; Vinyard et al., 2019). 

The regional differences observed in bonnetheads are 
greater than those that have been described for any 
adjoining populations of other elasmobranch species. Sig- 
nificant differences in growth characteristics of bonnet- 
heads have been found between populations in the GOM 


' Carlson, J. K., and J. Loefer. 2007. Life history parameters for 
Atlantic sharpnose sharks, Rhizoprionodon terraenovae, from the 
United States South Atlantic Ocean and northern Gulf of Mexico. 
Southeast Data, Assessment, and Review SEDAR13-DW-08, 7 p. 
[Available from website.] 


and the Atlantic region (Frazier et al., 2014). For the pop- 
ulation in the Atlantic region, estimated maximum age 
(males: 12.0 years; females: 17.9 years) and age at 50% 
maturity (males: 3.9 years; females: 6.7 years) are more 
than twice the estimates for the population in the GOM 
(males: 5.5+ years for maximum age, 1.7 years for age at 
50% maturity; females: 7.5+ years for maximum age, 2.9 
years for age at 50% maturity) (Lombardi-Carlson et al., 
2003; Frazier et al., 2014). Previous age and growth stud- 
ies of bonnetheads in the eastern GOM (Parsons, 1993; 
Carlson and Parsons, 1997; Lombardi-Carlson et al., 
2003) found significant latitudinal variation in life his- 
tory traits; however, latitudinal variation has not been 
detected in the population off the Atlantic coast (Frazier 
et al., 2014). High degrees of site fidelity have been doc- 
umented for this species (Heupel et al., 2006; Driggers 
et al., 2014), and tagging data indicate that there is no 
mixing between populations in the GOM and Atlantic 
region (Kohler and Turner, 2019). 

The population status of bonnetheads in U.S. waters 
was most recently assessed as a single stock in 2013; how- 
ever, because of observed differences in life history, tagging 
data, and genetic population structure, the results of that 
stock assessment were rejected and it was recommended 
that regional populations (i.e., those in the GOM and 
Atlantic region) be assessed as separate stocks (SEDAR, 
2013). As such, the population status for both stocks is cur- 
rently considered to be unknown (SEDAR, 2013). 

Recent studies of bonnethead population structure 
(Escatel-Luna et al., 2015; Portnoy et al., 2015) found 
that the populations in the GOM and Atlantic region are 
genetically distinct, with evidence of fine-scale genetic 
structure within populations. Using a combination of 
mtDNA and nuclear single nucleotide polymorphisms, 
Portnoy et al. (2015) found evidence of female philopatry 
with male-mediated gene flow. Results of further analysis 
indicate that over half of a small sample of outlier single 
nucleotide polymorphism loci has signatures of latitudinal 
selection. Portnoy et al. (2015) proposed that philopatry 
can lead to adaptive variation on a local scale and that, 
when combined with sex-biased dispersal, adaptive vari- 
ation can move among locations and environments. The 
high degree of site fidelity and latitudinal variation in life 
history observed for bonnetheads, in addition to localized 
adaption to environmental conditions, could explain the 
dissimilarities in life histories between populations in 
the different regions. However, differences in aging tech- 
niques, age estimation, or spatiotemporal biases could also 
explain observed differences (Campana, 2001; Cailliet and 
Goldman, 2004). 

The objectives of this study were 1) to use 2 long-term 
mark-and-recapture data sets and an age-independent 
model, GROTAG, to estimate region-specific growth rates; 
2) to generate estimates for age-independent life history 
parameters and compare results to region-specific esti- 
mates based on length-at-age data for verification of current 
life history information; and 3) to estimate region-specific 
seasonal growth and growth variability of bonnetheads in 
the GOM and the Atlantic region. 
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Materials and methods 
Data collection 


Data sets from long-term mark-and-recapture studies 
were available from 2 regions (the northeastern GOM and 
estuarine waters of the Atlantic coast of the southeastern 
United States), corresponding to the areas within each 
region where published, region-specific age and growth 
studies occurred. Mark-recapture data were provided for 2 
surveys conducted in the GOM. Mote Marine Laboratory 
provided data from fishery-independent gill-net surveys 
conducted from 1993 through 2006 in the eastern GOM 
(primarily from Yankeetown to Charlotte Harbor, Flor- 
ida). Detailed descriptions of the survey methods used 
by Mote Marine Laboratory can be found in Hueter and 
Tyminski’ and in Hueter and Tyminski (2007). Data for 
the northeastern GOM were obtained for the period from 
2003 through 2014 from the National Marine Fisheries 
Service from its Gulf of Mexico Shark Pupping and Nurs- 
ery (GULFSPAN) survey. A fishery-independent survey 
conducted with gill nets made of multiple stretched-mesh 
panels, the GULFSPAN survey is used to assess popula- 
tions of juvenile sharks from Cat Island, Mississippi, to 
Anclote Key, Florida, from April through October each 
year. Additional details about the GULFSPAN survey can 
be found in Bethea et al. (2015). 

For the population in the Atlantic region, bonnethead 
mark-recapture data collected from 1998 through 2019 
were available from the South Carolina Department of 
Natural Resources as part of its Cooperative Atlantic 
States Shark Pupping and Nursery (COASTSPAN) sur- 
vey. The COASTSPAN survey is conducted in estuarine 
waters along the coast of South Carolina from Saint 
Helena Sound to Bulls Bay during April-September each 
year. Detailed descriptions of COASTSPAN survey meth- 
ods can be found in Ulrich et al. (2007). 

For all surveys, upon capture, the precaudal length, fork 
length (FL), total length, and stretch total length of each 
shark were measured in a straight line along the axis of 
the body to the nearest half centimeter (for sharks from 
the GOM) or nearest millimeter (for sharks from the 
Atlantic region). If healthy, sharks were tagged externally 
either with a nylon dart tag (142-mm tag, Hallprint Fish 
Tags*, Hindmarsh Valley, Australia) or a FT-1-94 or T-bar 
Anchor Tag (Floy Tag Inc., Seattle, WA), at the base of the 
first dorsal fin (for both populations: GOM and Atlantic 
region), or with a 3.5-cm rototag (Dalton ID Systems Ltd., 
Henley-on-Thames, UK), inserted through the cartilage 
of the leading edge of the first dorsal fin (for the popula- 
tion in the Atlantic region only). A limited number of indi- 
viduals in the Atlantic region were double tagged with a 


” Hueter, R. E., and J. P. Tyminski. 2002. U.S. shark nursery 
research overview, Center for Shark Research, Mote Marine 
Laboratory 1991-2001. Mote Mar. Lab. Tech. Rep. 816, 31 p. 
[Available from website.] 
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rototag and a 12-mm, 125-kHz internal passive integrated 
transponder (Biomark Inc., Boise, ID) at the base of the 
first dorsal fin prior to release. 

Upon recapture during a survey, the date of capture, 
tag number, capture location (latitude and longitude), and 
all aforementioned length measurements were recorded 
before a shark was rereleased. If a bonnethead was recap- 
tured by a commercial or recreational fisherman (not 
during a survey), lengths, tag number, date of capture, and 
general capture location were requested. A subset of bon- 
netheads recaptured in the Atlantic region were sacrificed 
as part of age validation research (Frazier et al., 2014). 


Effects of tagging on growth 


Recaptured bonnetheads sacrificed and aged for age- 
based modeling of growth in a previous study (Frazier 
et al., 2014) were used in this study to test if there was 
any effect of tagging and tag type on growth. The residu- 
als at recapture (i.e., at sacrifice, calculated as expected 
FL at age minus observed length at estimated age) used 
for analysis were from use of the final VBGF in Frazier 
et al. (2014). Residuals at initial capture (i.e., at tagging) 
were determined by using FL at initial tagging and esti- 
mated age at initial tagging (calculated as estimated age 
at sacrifice minus time at liberty) to remodel the data set. 
The change in residuals at initial capture and in resid- 
uals at recapture of sacrificed bonnetheads were plotted 
against time at liberty. If tagging or tag type had a nega- 
tive effect on growth, most data points would be less than 
zero with the slope of the trendline significantly differ- 
ent from zero, indicating slower than predicted growth 
in tagged sharks. 


Modeling growth with tag-recapture data 


Growth increment data were modeled by using the Francis 
(1988b) method (i.e., by using the GROTAG model). In the 
event that a single individual was recaptured multiple 
times, only data from the initiai capture and final recap- 
ture were used in analyses to give equal weight to each fish 
and to maximize time at liberty (Welsford and Lyle, 2005). 
Data for all recaptured sharks (excluding data as previ- 
ously mentioned for sharks recaptured multiple times) 
were used in the model regardless of time at liberty, neg- 
ative growth, or potential outliers because the GROTAG 
model can use these data to inform several calculated 
parameters. Growth has been found to be significantly dif- 
ferent between sexes in previous growth studies (Parsons, 
1993; Carlson and Parsons, 1997; Lombardi-Carlson et al., 
2003; Frazier et al., 2014); therefore, sex-specific growth 
was modeled for both regions. 

The GROTAG model, which includes an implementa- 
tion of a maximum likelihood approach, was used to fit 
the VBGF (von Bertalanffy, 1938) to data for change in FL 
from initial capture to final recapture (AL) and for change 
in time at liberty (AT). The GROTAG model is a reparame- 
terization of the Fabens growth model (Fabens, 1965) that 
incorporates seasonal growth. Mean annual growth (g), 
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measured in millimeters FL per year, was estimated at 2 
user-defined reference lengths, « and B (where a<f). Ref- 
erence lengths were chosen to ensure that the majority of 
lengths at initial capture (L,) fell between the 2 defined ref- 
erence lengths by taking the mean length of the 3 smallest 
and 3 largest individuals at initial capture (a and B, respec- 
tively; Dureuil and Worm, 2015) in each region-specific data 
set. These parameters have better statistical properties 
than the asymptotic length, or the theoretical maximum 
length (L.,), and the growth constant (Rk) because they are 
not highly correlated. Further, they allow easier interpre- 
tation of growth from tagging data (Francis, 1988b). The 
growth rates relate to the parameters of the von Berta- 
lanffy growth curve as follows: 


8a. — & 
ek = 14 ot a SB (2) 
a—B 


Seasonal growth is parameterized as w, the time of year 
when growth is at its maximum, and as u, with a u value 
of 0.0 indicating no seasonal growth and a wu value of 1.0 
indicating strong seasonal growth with growth likely ceas- 
ing at some point during the year: 


8 AT +(¢>-6)) 
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where 6; = aa fori = 1, 2. 
Tt 


The GROTAG model is fit by minimizing the negative 
log-likelihood function (—A). Growth variability (v) is incor- 
porated into the model by the parameter wu; , the expected 
mean growth increment of the 7th individual where u; is 
normally distributed with a standard deviation (SD) of 
o;. In this study, 6, was assumed to be a function of the 
expected growth increment o,=vu;. An additional parame- 
ter p, the probability of outlier contamination, was also fit. 
For each data set, made up of i=1 to n growth increments 
where F is the range (largest and smallest) of observed 
growth increments, the following equation was used: 


Pp 
i= Byln|(.- ph +2] (4) 
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where i, = exp : 
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The likelihood function estimates the population measure- 
ment error in AL as being normally distributed with a mean 
of m and an SD of s. The initial model estimated g,,, 8, and v 
with additional parameters (m, s, w, uw, and p) added, increas- 
ing model complexity (Table 1). Unfitted parameters were 
held at zero. Optimal model parameterization was deter- 
mined for each region by using likelihood-ratio chi-square 
tests to determine if improvement in model fit was significant 
(P<0.05). Francis (1988b) suggested that the introduction of 


Table 1 


Parameters fitted for the GROTAG models used to esti- 
mate region-specific growth rates of bonnetheads (Sphyrna 
tiburo) in the Gulf of Mexico during 1993-2006 and in the 
Atlantic Ocean off the southeastern United States during 
1998-2019: growth rate estimates at reference lengths a 
and # (g, and gz), mean (m) and standard deviation (s) of 
the measurement error, magnitude (wu) and timing (w) of 
seasonal growth, growth variability (v), and outlier con- 


tamination probability (p). 


GROTAG model Parameters estimated 
Ew 8p, s 

Su &p> S,U 

Sa» &p, $V, mM 

S02 Ep» S, UV, M, U, W, P 
Ew &p> S,p 


an additional parameter should increase the log-likelihood 
value by at least 1.92. Likelihood ratio tests were also con- 
ducted to determine significant differences in von Bertalanffy 
growth curves between regions (Kimura, 1980). 


Age-based growth model 


To allow direct comparison of growth estimates based 
on age data and those based on tag-recapture data, 
region-specific length-at-age data from Lombardi-Carlson 
et al. (2003) and Frazier et al. (2014) were remodeled (the 
authors of both publications used the Beverton and Holt, 
1957, method of modeling VBGF parameters) by using an 
alternative parameterization of the VBGF recommended 
by Francis (1988a), in which mean length (L) of fish of age 
t is determined with this equation: 


(ly L lo “ 7 ve) 


L=l, + ae 


, (5) 


where p = 04x s 
ly —la 
where X = (® + Y)/2; 
l» = mean length at age ®; 
lx = mean length at age X; and 
ly = mean length at age V. 


Values for ® and were chosen to encompass the range 
of ages represented in the published length-at-age data 
from both regions (Lombardi-Carlson et al., 2003; Frazier 
et al., 2014). The Francis (1988a) parameterization yields 
estimates of VBGF parameters that better represent the 
growth information modeled from length-at-estimated- 
age and tagging data. The growth estimates generated 
from this model allow comparison of the mean growth rate 
of fish of an estimated age with that of fish of a length 
equal to the mean length at that age (Francis, 1988a). 
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Bootstrapping (with 5000 iterations) was used to develop 
95% confidence intervals (CIs) for VBGF parameter esti- 
mates. Model estimates and CIs were generated by using 
the FSA package (vers. 0.3.2; Ogle, 2012) in R (vers. 3.6.0; 
R Core Team, 2019). 


Comparisons of tag-recapture and age-based growth 


To compare growth rates between models, annual growth 
rates were estimated by using the parameters of the age 
and growth models to calculate mean length at estimated 
age. These values were used to create trend lines of the 
relationship between mean growth and FL, allowing 
comparison of the growth rates based on the GROTAG 
reference lengths. Parameter values from bootstrapping 
were used to calculate 95% CIs for the age-based growth 
rates. 

To allow statistical and visual comparison of growth 
rates between the populations in the Atlantic region and 
GOM, we calculated growth at 2 reference lengths corre- 
sponding to the mean of the 3 smallest L, of sharks from 
the Atlantic region and a length less than the smallest 
age-based bootstrap L., estimate (g,, gg). Region-specific 
bootstrap values (from the use of 5000 iterations) were 
used to calculate average growth and 95% CIs at the 
shared reference lengths. To allow comparisons of growth 
between age-based and length-based models, we calcu- 
lated mean growth rates for age-based models (g’,, g’3) 
corresponding to growth rates from length-based models 
(Sq. 8g). These were calculated by using bootstrapped esti- 
mates of parameters /», /x, and /y to solve for age-based 
growth rates (g’,, ’3). The resulting values equate to the 
expected annual growth rate of fish at the age when mean 
length is equal to reference length o or 8. The length-based 
(84, 8g) and age-based (g’,, g’3) growth rates were plot- 
ted with 95% CIs to allow comparisons of growth within 
and between regions and models. All GROTAG models 
were generated by using the fishmethods package (vers. 
1.4-0; Nelson, 2013) in R. 

For the GROTAG model, an additional parameter 
needed to produce von Bertalanffy growth curves, theoret- 
ical age at length zero (ty), was calculated by using mod- 
el-specific L., and k as well as observed mean size at birth 
(Lo) with age (t=0): 


1 (L,, — Ly) 
to = or Bu [a= (6) 


Region-specific mean Ly) was 263 mm FL for sharks from 
the Atlantic region (Frazier et al., 2014) and 252 mm FL 
for sharks from the GOM (R. Hueter, unpubl. data). 


Longevity and age at 50% maturity 


Estimated maximum ages from the region-specific age 
and growth studies provided initial values of maximum 
age; however, these values are likely underestimated in a 
fished population (Bishop et al., 2006). For the age-based 
and length-based models, longevity was estimated for 


each region and sex. The recaptured shark with the great- 
est time at liberty was assigned an age at initial capture 
by using L, and model-specific VBGF parameters recon- 
figured to solve for age where 


inf = 2 
L,, 


Age =| —— 


+ to. (7) 


This calculated age at initial capture was then added to 
time at liberty to estimate longevity. 

The GROTAG-derived VBGF parameters were used to 
solve for sex-specific age at 50% maturity by using equa- 
tion 7 with L, set to estimates of length at 50% maturity 
of sharks from the Atlantic Ocean (male: 617.8 mm FL; 
female: 818.5 mm FL; Frazier et al., 2014) and the GOM 
(male: 630.8 mm FL; female: 716.4 mm FL; Lombardi- 
Carlson‘), allowing calculation of age at 50% matu- 
rity without introducing bias from age _ estimation. 
Confidence intervals for longevity and age at 50% matu- 
rity were generated by using the 95% CIs for the VBGF 
parameters. For the published estimates of age at 
50% maturity, 95% Cls are recorded as published for the pop- 
ulation in the Atlantic region; however, no measurements of 
uncertainty were reported for the population in the GOM. 


Results 
Data collection 


A total of 139 and 190 recaptured bonnetheads were avail- 
able from the GOM and the Atlantic Ocean off the southeast- 
ern United States, respectively (Fig. 1). Ranges of lengths at 
capture and recapture and range and mean of time at liberty, 
by region and sex, are reported in Table 2. Sex- and region- 
specific recapture data sets did not encompass the range of 
lengths that have been reported as present in each region’s 
population (Ulrich et al., 2007; Bethea et al., 2015). Large 
females were lacking in the GOM data set (n=99, only 
2 females greater than 880 mm FL available; Suppl. Table 
[online only]), and there was no representation of small females 
in the Atlantic region data set (n=172, no females less than 
550 mm FL present; Table 2). Sample sizes for males from the 
GOM (n=40) and the Atlantic region (n=18) were relatively 
small, and lengths did not adequately represent the range of 
individuals found in the population, especially in the Atlantic 
region (Table 2, Suppl. Table [online only]). In comparison, the 
smallest and largest bonnetheads in the data used in the age- 
based models for fish in the Atlantic region are 245 and 
825 mm FL (males) and 262 and 1043 mm FL (females) 
(Frazier et al., 2014), and those in the data used in the GOM 


“ Lombardi-Carlson, L. A. 2007. Life history traits of bonneth- 
eads, Sphyrna tiburo, from the eastern Gulf of Mexico. South- 
east Data, Assessment, and Review SEDAR13-DW-24, 7 p. 
[Available from website.] 
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Figure 1 


Map showing the locations where bonnetheads (Sphyrna tiburo) were recap- 
tured in the northeastern Gulf of Mexico between 1993 and 2006 and in 
the estuarine waters of the Atlantic coast of the southeastern United States 
between 1998 and 2019. Circle size indicates the number of bonnetheads 
recaptured at a location. Data from these tag-recapture efforts were used in 
an age-independent model, GROTAG, to estimate growth rates of bonnet- 


heads by region. 


age-based models are 300 and 760 mm FL (males) and 340 
and 960 mm FL (females) (Lombardi-Carlson et al., 2003). 


Effects of tagging on growth 


A total of 22 recaptured individuals were aged by Frazier 
et al. (2014), and resulting data were available for analysis 
of the effect of tag type on growth (8 sharks were tagged 
with a nylon dart tag, and 14 sharks were tagged with a 
rototag). The changes in residuals from mean age at ini- 
tial capture (i.e., tagging) and estimated age at recapture 
were plotted against time at liberty (Fig. 2). The slope (6) 
of the line is not significantly different from zero for tag 
types combined or for individual tag types, and 95% CIs for 


Gullivan Baye Chokoloskee Bay 


slopes bound zero (all tags: b=0.005 [95% 
CI —0.033—0.043], P=0.615, df=20; nylon 
dart tag: b=0.037 [95% CI —0.109-0.184], 
P=0.554, df=6; rototag: b=0.005 [95% 
CI -0.042-0.052], P=0.805, df=12). There- 
fore, there is no evidence to indicate that 
tagging or tag type affected growth, and 
growth increments from tag-recapture 
data were considered suitable for model- 
ing growth in a population. 


Models based on tag-recapture data 


The GROTAG model (Francis, 1988b) 
produced biologically reasonable param- 
eter estimates for males from both 
regions (Table 3); however, 95% CIs are 
large, indicating that the sample size 
was insufficient to produce robust esti- 
mates of growth. The best-fit model for 
males in the GOM (model 3) included 
parameters for mean growth rates at 
reference lengths (g479, 8737), mean and 
SD of measurement error, and growth 
variability; the model failed to con- 
verge when the model 1 configuration 
was used (Table 4). The final model for 
males in the Atlantic region (model 2) 
included parameters for mean growth 
rates at reference lengths (8619, 765), 
SD of measurement error, and growth 
variability; model 3 failed to fit the data 
for this region (Table 4). Estimates of 
growth variability and seasonal growth 
for males are uninformative because the 
low sample size produced large 95% CIs 
for estimates, with upper and lower lim- 
its of 95% CIs above parameter bounds 
in the GOM model. Although confidence 
intervals are large and overlap, the SD of 
measurement error of the GOM model is 
3 times that of the Atlantic region model. 

Region-specific models for female bon- 
netheads converged for all models run. 
The final model for females in the GOM 
(model 2) included parameters for mean growth rates at 
reference lengths (g4¢5, 8915), 0D of measurement error, and 
growth variability (Table 4). The final model for females in 
the Atlantic region (model 4) is more complex, with the 
parameters for mean growth rates at reference lengths 
(555, £1000)) mean and SD of measurement error, growth 
variability, and seasonal variation included (Table 4). The 
estimates of growth variability for females are large from 
both the GOM model (v=0.63) and the Atlantic region 
model (v=0.56), indicating that individuals in the pop- 
ulation could be expected to grow 0.37—1.63 (GOM) or 
0.44-1.56 (Atlantic region) times the estimated average 
growth rate per length class. The model for females in the 
Atlantic region has a strong seasonal growth component 
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Table 2 


Summary of data collected for bonnetheads (Sphyrna tiburo) tagged and recaptured in the 
northeastern Gulf of Mexico (GOM) during 1993-2006 and in the Atlantic Ocean off the 
southeastern United States (Atlantic region) during 1998-2019. The number of individuals 
recaptured (7), the range of fork lengths (FLs) of fish at initial capture and recapture, and 
the range, mean, and standard deviation (given in parentheses) of times at liberty for recap- 


tured sharks are provided by region and sex. 


Mean time at 
liberty (d) 


Time at 
liberty (d) 


Initial Recapture 

Region n FL (mm) FL (mm) 
450-958 
520-780 
642-1014 
543-805 


1—2028 
1-1639 
10-3263 
13-2659 


GOM 99 430-930 
40 400-750 
Atlantic region 172 550-1013 
18 532-767 


308 (366) 
259 (319) 
458 (518) 
401 (823) 


(u=1.0, w=0.58), indicating that growth likely ceases at 
some point during the year and that it peaks in June. The 
SD of the measurement error for the model for females in 
the GOM (s=11.9 mm FL) is more than double that of the 
model for females in the Atlantic region (s=4.8 mm FL) 
with 95% CIs that do not overlap, indicating significant 
differences in measurement error between the 2 regions. 
The model for females in the Atlantic region had a negli- 
gible mean measurement error of —0.5 mm FL. The con- 
tamination probability parameter was not included in 3 
of the 4 models and is very low (p<0.001) in the model for 
females in the Atlantic region; therefore, the occurrence of 
outliers is scarce in all data sets. 


Age-based growth model and regional comparisons 


The Francis (1988a) age-based model produced nearly 
identical estimates of VBGF parameters to those produced 
by using the Beverton and Holt (1957) modeling method 
reported in Frazier et al. (2014). 

For both regions, the age-based models for males pre- 
dicted faster mean annual growth at smaller lengths and 
slower growth than the length-based models as bonneth- 
eads approached estimated L,, (Fig. 3). For males in the 
Atlantic region, 95% CIs for age-based estimates of growth 
rates overlap at GROTAG-predicted rates for both g¢1. 
and £75 (Fig. 3). For males in the GOM, 95% CIs overlap 
at GROTAG-predicted rates for 479; however, at 737, 95% 
CIs do not overlap and rates are significantly different 
from those from the age-based model, with the age-based 
model predicting near zero growth by an FL of 708 mm 
(Fig. 3). 

For females from the Atlantic region, the age-based 
model predicted a faster growth rate at smaller lengths 
but very similar growth rates at larger lengths compared 
with predicted rates from the GROTAG model (Fig. 4), 
and 95% CIs overlap at both g;,, and gj 999. For females 
in the GOM, the age-based model predicted nearly iden- 
tical growth rates at the smaller GROTAG reference 
length (g4¢5); however, the age-based model predicted a 


much smaller L,, than the estimate from the GROTAG 
model (Fig. 4). The predicted growth rate at the larger ref- 
erence length (go,;) is significantly higher than the esti- 
mate from the age-based model, with no overlap in 95% 
Cls. The significant differences in estimated growth rates 
for both males and females at g737 and go,;, between the 
GROTAG and age-based models indicate sampling bias or 
age underestimation in the GOM age-based model. 

Plots of growth rates and 95% CIs estimated for males 
with the GROTAG model do not indicate significant dif- 
ferences in growth between populations in the Atlantic 
region and the GOM (Fig. 3); however, results from like- 
lihood-ratio tests (x7=40.8, df=3, P<0.001) indicate signifi- 
cant differences in growth between regions. 

Plots of growth rates with 95% CIs estimated for 
females as well as results from likelihood-ratio tests 
(y7=31.2, df=3, P<0.001) for the best-fit GROTAG models 
indicate that growth was significantly different between 
regions, with a significant difference in average growth 
rates occurring in individuals larger than ~850 mm FL 
(Fig. 4). Plots of von Bertalanffy growth curves by sex, 
region, and model further illustrate the difference in pre- 
dicted lengths at age (Fig. 5); however, caution should be 
used when comparing curves between the GROTAG and 
age-based models because of different definitions of L 
(Francis, 1988a). 

To allow comparisons of growth within and between 
models and regions, growth rates were calculated for 2 ref- 
erence lengths (g;;,, 839) by using the GROTAG and age- 
based models. Plots of bootstrap parameter estimates 
from the GROTAG model indicate clear differences in g39 
between regions with no overlap in 95% Cls (Fig. 6). Plots 
of GROTAG bootstrap parameter estimates indicate less 
variation in growth with greater length for both regions. 
Overall, given the estimates from both models, the varia- 
tion in predicted growth is much higher for the population 
in the GOM than for the population in the Atlantic region. 
At the set reference lengths, growth rates do not signifi- 
cantly differ between the age-based and GROTAG models; 
however, as presented in Figure 3, estimated growth rates 
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Changes in residuals for ages of recaptured bonnetheads (Sphyrna tiburo) in 
the estuarine waters of the Atlantic coast of the southeastern United States 
plotted against their times at liberty, by type of tag: nylon dart (squares) and 


3000 3500 


y=0.0054x— 14.886, R?=0.0053 


overlap (Table 3). For sharks in the GOM, 
published estimated maximum ages are 
significantly below calculated maximum 
ages from models based on age and on 
tag-recapture data, with published val- 
ues below 95% Cls for all models with the 
exception of the GROTAG model used for 
estimating growth of males; a lower limit 
of the 95% CI could not be estimated 
because of the size of the tagged shark 
used in relation to the lower limit of the 
95% CI for the estimate of L,,. For bonnet- 
heads in the Atlantic region, calculated 
maximum ages are also significantly 
larger than published estimates, with the 
lower limit of the 95% CI falling above 
maximum estimated age for all models. 

Published estimates of age at 50% 
maturity are lower than those calculated 
by using the GROTAG models (Table 3). 
We were unable to calculate the lower 
limit of the 95% CI for age at 50% matu- 
rity of males with either GROTAG 
model, and 95% CIs based on the GOM 
data were not published; therefore, we 
cannot determine if differences are sig- 
nificant for these estimates. Upper and 
lower confidence intervals for estimates 
from both length-based and age-based 
models were available only for females 
in the Atlantic region. Although the 
GROTAG model for females in the Atlan- 
tic region produced an estimate of age at 
50% maturity that is slightly larger than 
that from the age-based model, 95% CIs 
for the estimates overlap. 


rototag (triangles). If tagging has a negative effect on growth, the majority 
of data points should be negative and the slope of the dashed line should be 


significantly different from zero. The equation for the trend line and the coef- 


Discussion 


ficient of multiple determination (R) are provided. Residuals were calculated 


by using lengths at initial tagging and at recapture and by using ages and the 


von Bertalanffy growth curve from Frazier et al. (2014). 


deviate between models as individuals grow beyond the 
upper reference length. 


Longevity and age at 50% maturity 


For the population in the GOM, data for recaptured bonnet- 
heads from long-term tag-recapture studies, a male and a 
female with times at liberty of 1639 and 1833 d and initial 
FLs of 600 and 800 mm, respectively, were used to calcu- 
late maximum age in each model. For the population in the 
Atlantic region, data for a recaptured male and a recap- 
tured female, with times at liberty of 2659 and 3263 d and 
FLs of 763 and 980 mm, respectively, were used. Calculated 
ages from the GROTAG model are larger than those from 
the age-based models for both regions, although 95% CIs 


This study confirms previously pub- 
lished significant differences in region- 
specific growth and life history 
characteristics of bonnetheads between 
the northeastern GOM and the Atlantic region by using 
age-independent methods of modeling growth. Growth 
rates, age at 50% maturity, and longevity all differ for both 
sexes between these 2 regions. Despite low sample sizes and 
poor coverage of ranges of lengths for males, we were able 
to produce reasonable growth estimates for the populations 
in the northeastern GOM and the Atlantic region by using 
the GROTAG model; however, we were unable to determine 
significance of some life history characteristics because of 
the high uncertainty in estimates of growth and growth 
parameters. The low sample sizes of males from the Atlan- 
tic region available for growth modeling is reflective of the 
nearshore distribution of male bonnetheads, with tagging 
effort heavily skewed toward females in estuarine waters 
where most fishery-independent sampling and tagging 
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Table 3 


Estimates of von Bertalanffy growth function (VBGF) parameters, maximum age, and age at 50% maturity from GROTAG models 
based on tag-recapture data and from age-based models for male and female bonnetheads (Sphyrna tiburo) from the northeast- 
ern Gulf of Mexico (GOM) and the Atlantic Ocean off the southeastern United States (Atlantic region). The VBGF parameters, 
asymptotic length (L..), coefficient of growth (k), and theoretical age at length zero (¢)), are provided with 95% confidence intervals, 
which were generated by using bootstrapping (with 5000 iterations). Also provided are published estimates of maximum age 
(Lombardi-Carlson et al., 2007; Frazier et al., 2014). Ages at 50% maturity were estimated in GROTAG models by using lengths at 
50% maturity published in Lombardi-Carlson et al. (2007) and Frazier et al. (2014). Data used in the models are from sharks tagged 
and recaptured in the GOM during 1993-2006 and in the Atlantic region during 1998-2019. FL=fork length. 


Model 


GROTAG 


Atlantic region 


Age based 


Atlantic region 


L., 
(mm FL) 


769.6 
604.2—823.5 
948.3 
921.2—967.9 
814.8 
533.4-1135.6 
1039.5 
1027.7—1048.5 
703.3 
663.8—765.2 
895.0 
842.9-970.3 
7179.8 
761.5—799.7 
1032.3 
1011.6—1053.9 


k 


0.254 
0.470-0.087 
0.243 
0.313-0.184 
0.166 
0.386—0.070 
0.170 
0.191-0.145 
0.538 
0.753-0.386 
0.282 
0.334—0.226 
0.296 
0.320-0.274 
0.188 
0.179-0.198 


to 
(years) 


-1.56 
—(1.14—4.19) 
-1.27 
—(1.02-1.64) 
—2.35 
—(1.76-3.76) 
-1.75 
—(1.54-1.99) 
-1.60 
—(1.10—2.25) 
-2.13 
—(1.82—2.55) 
-1.51 
—(1.35—-1.68) 
-1.76 
—(1.59-1.94) 


Age at 50% 


maturity 
(years) 


5.2 
oo12.5 
4.5 
3.85.7 
6.2 
oo—7.5 
7.4 
6.8-8.5 
2.0+ 


3.0+ 
3.9 


6.7 


Maximum 


age 
(years) 


8.9 
oo—15.3 
11.4 
9.7-14.3 
21.5 
19.4—co 
24.0 
23.5-25.8 
6.5 
6.4-6.6 
10.8 
10.2-12.1 
18.7 
16.8—co 
23.0 
21.9-24.9 


Published 
maximum 


age 
(years) 


occurs (Ulrich et al., 2007). In the GOM, the lower num- 
ber of male bonnetheads that were recaptured is a result 
of the primary habitat of male bonnetheads being outside 
of or on the periphery of the area where sampling occurred 
and where tags were deployed. Therefore, fewer males than 
females were tagged and consequently recaptured. Future 
efforts should focus on sampling and tagging male bonneth- 
eads, especially those in early life stages, to decrease uncer- 
tainty in results from length-based models. 

Although it is difficult to determine the effects of tagging 
on growth in wild fish populations, our results do not pro- 
vide evidence indicating that tagging and tag type affected 
growth. Growth was variable in the first year after tagging, 
with a larger number of individuals experiencing slower 
growth than expected. This slow growth, however, could be 
due to the short-term effects of the stress from being cap- 
tured rather than the effect of the tag on the individual 
(Gruber, 1982; Parsons, 1987; Davenport and Stevens, 1988; 
Skomal and Bernal, 2010). Over the long term, we found no 
evidence of effects of tagging on growth of captured bonnet- 
heads, and all individuals at liberty for over 4 years grew 
older and larger than the age and length predicted by the 
age-based growth model. These findings contrast with those 
of studies on lemon sharks (Negaprion brevirostris) (Manire 
and Gruber, 1991; Oliveira, 2001) and northern pike (Esox 
lucius) (Scheirer and Coble, 1991) that indicate that growth 


could be impaired by up to 50% in tagged individuals when 
a variety of tag types were used. Results of other growth 
studies (e.g., Jensen, 1967; Jolivet et al., 2009) indicate 
that there were no significant effects of tagging on growth; 
therefore, it is possible that there may be effects specific 
to species, life stage, or tag type. Caution should be used 
in interpreting growth data from recaptured individuals to 
the population level without investigating these effects. 
The final GROTAG models for males in both regions 
included SD of measurement error and growth variability, 
and the GOM model has the added parameter of mean 
measurement error. Estimates of growth variability are 
unreliable because of low sample size, and 95% CIs con- 
tain values that are not above and below parameter limits; 
therefore, no comparisons can be made with these param- 
eters. The final GROTAG models for males and females in 
the GOM both have large estimates of SD of measurement 
error compared with those from the models for males and 
females in the Atlantic region; however, this result was not 
unexpected. For individuals from the GOM, length was 
measured to the nearest half centimeter, and several biol- 
ogists and interns measured fish; whereas, bonnetheads 
from the Atlantic region were measured to the nearest 
millimeter by only 2 individuals. The use of fewer individ- 
uals to measure all captured specimens and the smaller 
measurement increment likely led to more precise 


Fishery Bulletin 118(4) 


338 


VST 0 
9¢6T 0 


98€'0-040'0 
99T'0 
€LT 0 
S810 
S810 
0LV'0—-L80'0 
¥Gc 0 
9960 


69T'0 


9°S68 

€'668 

9 SETI-V' EES 
871s 

6618 

COLL 

COLL 
G'E68-ZG F09 
9'°69L 

9'G8L 


V 6801 


I6L0-G7T0 G8POI-L LCOT 


LOT'O 
OsT0 
TLT0O 
0660 
€9T0 
866 0 
€66'0 
€1€0-78T 0 
676 0 
9760 


4 


G'6E0T 

9 EPOT 
L801 

6 O€0T 
8°G96 

G6 1S6 

T 096 

6 196-6 166 
€ 876 

0876 


(TH wut) 
“1 


sreqyoueied TOCA 


GTVG- 
g'0- 
6G 


= OS 

0:0 eV 

0S'0-900- Ga8-ZGY 

GG0 LG 

= G8 

= € 61 

OT 5 A0) gv 
= 08 T-00T VPE-O'L 
= 00'T 6ST 
= 00'T 9ST 


= = = U8 
g9°0-970 OO0T-780 90-970 OL-6Z 
890 OT 99°0 8v 
= a 09°0 L8 
TS'0 86 

= GVG 

a GLI 

v9'0 9TT 

L9°0 STI 

180-970 TST-8'8 

€9'0 6 TT 

= GCE 


(‘Ta uu) 
s 


siloyouered pajepoyy 


L’8 9°0€ 

OL 8'SG 
86-TG- €8P-0CT 
91 6 0€ 

OL 6 TE 

9g 10S 

9G L'0S 
GI¢-6y- 8YII-OFT 
EL €L9 

Or LOL 

T9 6 PL 
€L-G7 L’€8-T'99 
T9 LUL 

T9 €'89 

T9 GOL 

T9 666 

GL 67L 

GL 0°66 

GL v'66 
VOI-€&E Y98zGI-G'ES 


GL € VOL 
L'6 Vv SOT 


dg 


(aead/T Aq UU) 
SO7CUITISO YMOLD 


pooyrpeyty 


WOD 9[eWey 
[epoy uote lal 
(Ta wu) 


-30T yysue] 
aoUdIIJOY 


YASUE] YLQJ="TI ‘900S-E66T SuLnp uvsoQ onUELIY 


oY} UL pUS 900Z-E66T SULINP JOD 0} Ul poinydeoer pue passe syeYs WOIJ ole S[aPOU dy UI pasn eyeC ‘(,.) SYSLI9}Se YYIM poyeoIpUr are YoY ‘sjepour [eu 10; pozzodes 
a1 S[BAIOJUI BDUEpPYUOI % GG oy, ‘poyzoder ore (Y) JUSTOTJJOOO YYMOId pue (~'7) YASUE] o1y03duIASe Jo Si1ojoweIed (qO_A) UorjUN; YyMOIS AjyULTeJog UOA pue ‘(d) AyIqeqord 
UOTJEUIUTLJUOD IaTTINO “(a) APTIGVIIeA YYMOIS ‘YZMOIS [eUOSveS Jo (2) ZUTUIT] pue (7) OpnyzUseUI “1OLI JUSULOINSvOUT OY} Jo (S) UOIYeIASp paepueys pue (uw) uvou ‘(IF pue 9) 
SYISUS| SUAIEJOL eV So}eI YYMOIS JO SoyeUIT}sSa ‘(Y) ONTeA pooyl[exl]-so] ‘(¢] pue 0) Japoul Yove AO} SYASUI] sdUDIBJoy “UOTSeI pue xas Aq ‘(‘[}V) S9}eIS popup UTEYseoy4NOs oy} 
JJO UBIO ITJULTY OY} UI pue (OY) OdTxeJ[ JO JIND U1S}SeoyyAOU 9Y] UT YYSNVd (0uNG?? DULMYdgG) speoyjouUog Oj eYep o1n4dedoI1-3e} UO pase S[Epoul HY LOWD Woy szpnsey 


v P1921 


Frazier et al.: Growth rates of Sphyna tiburo estimated from tag-recapture data 


-e— GROTAG model 
-& Age-based model 


100 120 
nl 


40 60 80 
L 


20 


-e- GROTAG model 
-& Age-based model 


80 100 120 


60 


40 


Growth rate (mm FL/year) 
20 


-4- Atlantic region 
-~ GOM 


So 
=) 
i 
g 
=) 


60 


20 40 


0 


550 650 
Fork length (mm) 


Figure 3 


Growth rates of male bonnetheads (Sphyrna tiburo) estimated with the opti- 
mal GROTAG models, which are based on tag-recapture data, and mean 
growth rates calculated with the age-based (Francis, 1988a) von Bertalanffy 
growth function, with associated 95% confidence intervals (dotted lines). 
Growth rates are shown for bonnetheads tagged (A) in the northeastern Gulf 
of Mexico (GOM) and (B) in the estuarine waters of the Atlantic coast of the 
southeastern United States (Atlantic region). Also provided are (C) estimates 
of growth rates from the optimal GROTAG model for each of the 2 regions. 
Age-based models were generated by using original length-at-age data from 
Lombardi et al. (2007) for sharks in the GOM and from Frazier et al. (2014) 
for sharks in the Atlantic region. Data used in the GROTAG models are from 
bonnetheads tagged and recaptured in the GOM during 1993-2006 and in the 
Atlantic region during 1998-2019. Error bars indicate 95% confidence inter- 
vals generated by bootstrapping (with 5000 iterations). FL=fork length. 
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and differences in these growth rates 
may be an artifact of sample size given 
the observed 95% CIs. Regardless, these 
variability estimates are among the 
highest published values for elasmo- 
branchs with values ranging from 0.06 
for the tope (Galeorhinus galeus) (Fran- 
cis and Mulligan, 1998) to 0.58 for the 
spotted estuary smooth-hound (Mustelus 
lenticulatus) (Francis and Francis, 1992), 
with most falling in a range of 0.15—0.40 
(e.g., Simpfendorfer, 2000; Meyer et al., 
2014; Natanson and Deacy, 2019). 
Results from the final GROTAG mod- 
els indicate a strong seasonal growth 
pattern in the population in the Atlantic 
region (u=1) but not in the population in 
the GOM (the seasonal growth parame- 
ter was not included in the final GOM 
model). Latitudinal variation in growth 
has been previously detected for bon- 
netheads in the GOM but not for those 
in the Atlantic region (Parsons, 1993; 
Carlson and Parsons, 1997; Lombardi- 
Carlson et al., 2003; Frazier et al., 2014). 
The more variable growth for sharks 
in the GOM and lack of a seasonal sig- 
nal may indicate that there are differ- 
ences within the population in growth 
along the western coast of Florida, 
where sampling for this study occurred. 
Site fidelity has been established for 
bonnetheads in both the Atlantic region 
and the GOM (Heupel et al., 2006; 
Driggers et al., 2014), and bonnetheads 
are known to make predictable tem- 
perature-driven migrations in response 
to declining water temperatures (Ulrich 
et al., 2007; Driggers et al., 2014). In the 
Atlantic region, all bonnetheads were 
tagged and recaptured within their 
established summer range. In the GOM, 
bonnetheads were also predominantly 
captured and recaptured during sum- 
mer months; however, capture and 
recapture of sharks took place at multi- 
ple locations along the coast. Because 
these sampled subpopulations likely 
have site fidelity to their area of cap- 
ture, they likely experience differential 
growth due to differences in water tem- 
peratures (Carlson and Parsons, 1999; 
Pistevos et al., 2015), in food availabil- 


measurements at tagging and recapture and consequently 
lower chances of errors in measurement. 

Results from the final GROTAG models for females indi- 
cate that growth is more variable for bonnetheads from 
the GOM (v=0.63) than for those from the Atlantic region 
(v=0.56); however, these differences are not significant 


ity or energetic value of food (Bethea et al., 2007; 
Vucic-Pestic et al., 2011), and in energetic needs during 
winter migrations necessary to find optimum water tem- 
peratures for overwintering (Carlson and Parsons, 1999; 
Hoffmayer et al., 2006). However, the lack of a detection 
of seasonal growth patterns and greater growth 
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Figure 4 


Growth rates of female bonnetheads (Sphyrna tiburo) estimated with the 
optimal GROTAG models, which are based on tag-recapture data, and mean 
growth rates calculated with the age-based (Francis, 1988a) von Bertalanffy 
growth function, with associated 95% confidence intervals (dotted lines). 
Growth rates are shown for bonnetheads tagged (A) in the northeastern Gulf 
of Mexico (GOM) and (B) in the estuarine waters of the Atlantic coast of the 
southeastern United States (Atlantic region). Also provided are (C) estimates 
of growth rates from the optimal GROTAG model for each of the 2 regions. 
Age-based models were generated by using original length-at-age data 
from Lombardi et al. (2007) for sharks in the GOM and from Frazier et al. 
(2014) for sharks in the Atlantic region. Data used in the GROTAG mod- 
els are from bonnetheads tagged and recaptured in the GOM during 1993-— 
2006 and in the Atlantic region during 1998-2019. Error bars indicate 95% 
confidence intervals generated by bootstrapping (with 5000 iterations). 
FL=fork length. 
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Results from GROTAG models con- 
firm previously published differences 
in growth between the 2 regions. Both 
males and females in the GOM had 
faster average growth at smaller lengths 
and slower growth with a smaller L. 
than bonnetheads from the Atlantic 
region. Visual inspection of 95% CIs 
does not reveal significant differences 
between regions because of large uncer- 
tainty caused by low sample sizes, but 
results of likelihood-ratio tests confirm 
differences in growth. Significant differ- 
ences in growth have been detected for 
other coastal migratory sharks (Carlson 
and Baremore, 2003; Driggers et al., 
2004; Vinyard et al., 2019); however, the 
reasons for these regional growth differ- 
ences are unknown. 

Growth curves derived from GROTAG 
and age-based models are not directly 
comparable; instead, Francis (1988a) 
suggested that comparisons of growth 
rates between data types are more 
appropriate. All GROTAG models pre- 
dicted slower initial growth and faster 
growth at longer lengths compared with 
estimates from age-based models. The 
slower estimated initial growth rates 
likely result from a lack of tagged young- 
of-the-year bonnetheads in all 4 data 
sets, likely a combination of gear bias 
and variation in spatial distribution of 
young-of-the-year bonnetheads (Drig- 
gers et al., 2014). The GROTAG model for 
females in the GOM is based on the larg- 
est sample size of tagged and recaptured 
juvenile bonnetheads (including indi- 
viduals tagged as small as 400 mm FL) 
and, consequently, the best agreement 
between model growth rates at small 
lengths. However, growth rates are 
significantly different at long lengths. 
Estimates from the GROTAG model for 
females in the Atlantic region have excel- 
lent agreement with those from the age- 
based model, especially as lengths reach 
predicted L... Parameters of the VBGF 
are not significantly different between 
models, confirming estimated growth 
from the age-based model for females in 
the Atlantic region. Although the VBGF 
parameters are defined differently for 


the length-at-age data and tag-recapture data, the results 
of other studies indicate that models based on tag-recap- 
ture data can produce growth curves that are similar to 
those from age-based models (Natanson et al., 2002). 

The GROTAG models for males and females predicted 
significantly faster growth at longer lengths, larger L.,, 


variation may also be artifacts of sample size. Future 
efforts should focus on continued tagging of bonnetheads 
in these regions because a more robust data set would 
allow age-independent comparisons of growth and, there- 
fore, determination of latitudinal differences in growth 
within the GOM. 
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Figure 5 


Comparison of growth curves from GROTAG models, which are based on tag- 
recapture data, and from age-based von Bertalanffy growth functions based 
on data from Lombardi et al. (2007) and Frazier et al. (2014) for (A) female 
and (B) male bonnetheads (Sphyrna tiburo) tagged in the northeastern Gulf of 
Mexico (GOM) during 1993-2006 and off the Atlantic coast of the southeastern 
United States (Atlantic) during 1998-2019. Age-based models were generated 
by using original length-at-age data from Lombardi et al. (2007) for sharks in 
the GOM because published models used total length and band count. 
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model for females in the GOM may be an 
artifact of sample size, although the 
model produced much greater variation 
in growth estimates at gg3), contrary to 
the estimates from other models indicat- 
ing a source of variation other than sam- 
ple size. These results indicate that there 
may be significant underestimation of 
ages of bonnetheads in the GOM with 
the age-based models. 

Sampling, fishing pressure, or gear 
bias could have affected results of the 
age-based models; however, estimates of 
age at 50% maturity calculated by using 
age-independent, length-based growth 
estimates also point to significant dif- 
ferences in maturity. These estimates of 
age at 50% maturity indicate that female 
bonnetheads mature at an average of 
1.5 years and males mature at an average 
of 3.2 years later than females and males 
do according to age-based estimates, 
providing further evidence of age under- 
estimation in the GOM data set. The GRO- 
TAG model does not produce comparable 
estimates of average age at length; how- 
ever, data from long-term tag-recapture 
studies were available. When these data 
were combined with results from growth 
models, calculated longevity for female 
and male bonnetheads in the GOM are 
8.9 and 11.4 years from the GROTAG 
models and 6.5 and 10.8 years from the 
age-based models; these ages are con- 
siderably older than the maximum age 
estimates of 5.5+ and 7.5+ reported by 
Lombardi-Carlson et al. (2003). 

Data from long-term tag-recapture 
efforts and longevity estimates calcu- 
lated from GROTAG models confirm the 
report of age underestimation for bon- 
netheads in the Atlantic region by Fra- 
zier et al. (2014). However, for females, 
length-based estimates of age at 50% 
maturity are not significantly different 
from age-based estimates. If age under- 
estimation was common in the age-based 
study for sharks in the Atlantic region, it 


and faster rates at which k is reached compared with esti- 
mates from age-based models. Plots of bootstrapped 
growth estimates for females from the GOM age-based 
model indicate large uncertainty in growth estimates at 
the upper reference length. Comparatively, results from 
GROTAG models for bonnetheads in the Atlantic region 
and in the GOM and from the age-based model for bonnet- 
heads in the Atlantic region indicate less uncertainty at 
the upper reference length compared with the lower refer- 
ence length at which growth is expected to be more vari- 
able (Erzini, 1994). The large uncertainty in the age-based 


may be rare or only in individuals at or near L,, given the 
agreement between growth models and life history param- 
eters for females. 

Although age and growth studies are not necessarily 
expected to encounter the oldest individuals in a popula- 
tion because of their relative scarcity (Bishop et al., 2006), 
there is increasing evidence that age is frequently under- 
estimated in many studies of age and growth of elasmo- 
branch species, especially in those of long-lived species 
(e.g., Kalish and Johnston, 2001; Francis et al., 2007; 
Frazier et al., 2015; Harry, 2018; Natanson et al., 2018b). 
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more reliable, estimates of longevity 
and age-independent estimates of matu- 
rity. Although methods of modeling 
growth with tag-recapture data have 
historically been considered inferior to 
age-based methods (Harry, 2018), given 
recent doubts about the suitability of 
vertebrae for aging, methods based on 
tag-recapture data may be more reliable 
than age-based methods (Harry, 2018; 
Natanson and Deacy, 2019). When pos- 
sible, more computationally advanced 
methods, such as those of Aires-da-Silva 
et al. (2015) and Francis et al. (2016), 
which incorporate all available growth 
data (length based, tag and recapture, 
and age based) that can be integrated 
directly into stock assessment models, 
should be considered because they can 
combine all available data sources to 
estimate growth characteristics of a 
population. Additionally, methods that 
involve tag-recapture data may be useful 
in providing growth information for spe- 
cies that are difficult to age and for pop- 
ulations for which sacrifice of specimens 
is not preferred or fishery-dependent 
catch are not available for sampling (e.g., 
for protected species). 

Given the large discrepancies between 
the age-based model and the model based 
on tag-recapture data for bonnetheads in 
the GOM, the length-based tag-recapture 
estimates generated in our study (espe- 
cially age at 50% maturity and longevity) 
likely more accurately describe growth of 


Plot of bootstrap estimates of growth rates (g) at the 2 reference lengths of 
555 and 830 mm fork length (FL), parameters g;,, and 83), from GROTAG 
models based on tag-recapture data for female bonnetheads (Sphyrna tiburo) 
tagged in the northeastern Gulf of Mexico (GOM) during 1993-2006 and off 
the Atlantic coast of the southeastern United States (Atlantic region) during 
1998-2019 and from region-specific, age-based von Bertalanffy growth func- 
tions generated by using original length-at-age data from Lombardi et al. 
(2007) for sharks in the GOM and from Frazier et al. (2014) for sharks in the 
Atlantic region. Mean growth rates for age-based models were calculated from 
bootstrapped estimates of lengths at ages ®, X, and Y (parameters lg, ly, ly). 
Estimates are shown with 50%, 80%, and 95% polygon ellipsoid confidence 
intervals. Bootstrapping was done with 5000 iterations. 


bonnetheads in the GOM than age-based 
estimates. Parameters of the VBGF are 
commonly used in modeling estimates 
of mortality, maturity, and longevity 
and in modeling for stock assessments. 
Use of current age-based parameters 
for the population in the GOM should 
be carefully considered until growth in 
the population can be reexamined. Using 
parameters that underestimate age and 
maturity can lead to overly optimistic 
estimates of growth and mortality rates, 
increasing the potential for overexploita- 
tion of stocks. Results of this study indi- 


Results of recent studies indicate that in some species of 
sharks, band-pair formation may be a function of growth 
and body girth, not age (Natanson et al., 2018b). The 
results of our study indicate that the use of tag-recapture 
methods can produce estimates of growth that are similar 
to those produced with age-based methods. 

When sufficient tag-recapture data are available, mod- 
els should be used to verify growth estimates from age- 
based models and provide independent, and perhaps 


cate that bonnetheads in the GOM may be more susceptible 
to overexploitation than previously believed. The status of 
the stocks of bonnetheads in the northeastern GOM and 
off the Atlantic coast of the southeastern United States are 
currently considered to be unknown, and no stock assess- 
ments are pending. Given these results, we recommend 
that region-specific assessments be conducted and that 
they incorporate these new age-independent estimates of 
life history characteristics. 
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Abstract—Effective assessment and 
management of migratory species 
require an understanding of stock 
boundaries and habitat use during crit- 
ical life history periods. The blue crab 
(Callinectes sapidus) is considered a 
primarily estuarine species and has 
traditionally been assessed and man- 
aged on a state basis, relying on abun- 
dance estimates derived from estuarine 
sampling. Females undertake a sea- 
ward spawning migration, and recent 
evidence indicates that coastal and 
offshore waters may provide import- 
ant spawning habitat. We investigated 
occurrence and reproductive output of 
spawning blue crab in offshore waters 
of the northwestern Gulf of Mexico, 
and results indicate that these waters 
may be an important spawning area 
for this species. Densities of blue crab 
were low offshore compared with those 
in estuaries, but given the large extent 
of available offshore habitat, the total 
abundance of adult females is likely 
similar in offshore and estuarine 
waters. Analyses of the reproductive 
output of females collected offshore 
indicate that they continue to spawn 
after leaving the estuaries. Fecundity 
and measures of egg quality were 
assessed and were similar to previous 
estimates for estuarine regions. Future 
management and assessment efforts 
should account for this previously 
understudied portion of the blue crab 
spawning stock. 
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Fisheries management is most effective 
when the scale of management matches 
the scale of the stock (Begg et al., 1999; 
Al-Humaidhi et al., 2013). The modern 
stock concept describes units of a pop- 
ulation that can be considered homoge- 
neous (in terms of life history, responses 
to environmental fluctuations, genetic 
diversity, etc.) and can inform the scale 
of assessments and management deci- 
sions (Begg and Waldman, 1999). The 
failure to correctly identify spatial 
complexity within a stock has been 
identified as one of the main factors 
contributing to the collapse of some 
high-profile fisheries in the recent 
past (Hutchinson, 2008), and studies 
examining genetic stock structure have 
identified cases with clear mismatches 
between the genetic stock structure 
and the management of the stock as a 
singular unit (Bonanomi et al., 2015; 
Kerr et al., 2017). 

Migratory stages are common in 
marine and estuarine species (e.g., 
Forward and Tankersley, 2001; Potter 
et al., 2015; Secor, 2015). Large-scale 
seasonal or annual movements can 
complicate fisheries management 
because different life stages may 
occupy different habitats. For these 


species, understanding habitat use 
during critical life history stages (e.g., 
spawning) is particularly important 
in the context of temporal or spatial 
management actions (e.g., season and 
area closures, spawning sanctuar- 
ies, migratory corridors, and ecosys- 
tem-based management). 

Blue crab (Callinectes sapidus) sup- 
port highly valuable fisheries along the 
Atlantic and Gulf coasts of the United 
States, with landings in 2018 totaling 
over 62,000 metric tons and having a 
wholesale value of over $193 million 
(NMFS!). Following their terminal 
molt, female blue crab mate and store 
the sperm necessary for their life-long 
reproductive output (Millikin and 
Williams, 1984). After mating, female 
blue crab forage for several weeks 
(Turner et al., 2003) before begin- 
ning a seaward spawning migration 
(Aguilar et al., 2005). This migration 
ensures that females reach spawning 
grounds of sufficiently high salinity 
(>20.1; Costlow and Bookhout, 1959) 


' NMFS (National Marine Fisheries Service). 
2020. Fisheries of the United States, 2018. 
NOAA, Natl. Mar. Fish. Serv., Curr. Fish. 
Stat. 2018, 140 p. [Available from website.] 
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for larval survival. Spawning female blue crab are often 
found in large numbers at or near the mouths of estuaries 
(e.g., Rittschof et al., 2010; Anderson et al., 2017); these 
areas are typically assumed to be the primary spawning 
grounds. The blue crab is a highly fecund species, with 
females producing multiple clutches of several million lar- 
vae at a time (Hines et al., 2003). Each brood of larvae 
is released into the water column as zoeae that disperse 
through currents for 30-50 d (Costlow and Bookhout, 
1959) before their return and recruitment to coastal habi- 
tat as megalopae soon after (Rabalais et al., 1995; Ogburn 
et al., 2009). 

Mature females can travel great distances (>200 km; 
Aguilar et al., 2005) during their spawning migration, 
and migration rates can exceed 8 km/d (Carr et al., 2004). 
Female blue crab continue migrating throughout their 
reproductive period (Hench et al., 2004) and are, there- 
fore, likely moving out of estuaries and continuing their 
migration in offshore waters. Adults of this species typ- 
ically are considered estuarine, but large numbers of 
spawning female blue crab have been found in waters off 
both the Atlantic coast (Rittschof et al., 2010; Ogburn and 
Habegger, 2015) and Gulf coast (Gelpi et al., 2009, 2013) 
of the United States. Although observed densities of blue 
crab are typically lower offshore than in estuaries, evi- 
dence from the Atlantic coast of the southeastern United 
States indicates that this often-overlooked offshore popu- 
lation may represent a substantial fraction of the regional 
spawning stock (Ogburn and Habegger, 2015). In the 
Gulf of Mexico, Gelpi et al. (2009, 2013) found large num- 
bers of spawning females residing on Ship, Trinity, and 
Tiger Shoals and surrounding waters, 20-40 km off the 
coast of Louisiana. All developmental stages of eggs were 
observed, and 81% of females had full ovaries at the time 
of capture, indicating that they would continue to spawn 
at least one additional clutch of eggs during the season 
(Gelpi et al., 2009). Furthermore, spawning female blue 
crab are regularly captured in trawl surveys conducted 
by the Southeast Area Monitoring and Assessment Pro- 
gram (SEAMAP) in the Gulf of Mexico, as far as 200 km 
offshore. 

Given the large extent of relatively shallow offshore 
habitat in the Gulf of Mexico and with offshore habitat 
proving to be of larger importance to the spawning stock 
of blue crab in the Atlantic Ocean than previously thought 
(Ogburn and Habegger, 2015), we sought to identify the 
contribution of this offshore habitat to the spawning stock 
of blue crab in the northwestern Gulf of Mexico. The goal 
of this study was to conduct an investigation of spawning 
female blue crab in offshore waters of the northwestern 
Gulf of Mexico. We examined temporal and spatial varia- 
tion in catch of blue crab in fishery-independent SEAMAP 
groundfish trawl surveys over a 20-year period (2000- 
2019) and examined the reproductive output and overall 
condition of females collected in offshore waters during 
the SEAMAP groundfish trawl surveys in 2017. We also 
compared results for spawning females found offshore 
in our study with previous findings for blue crab found 
Spawning in inshore areas. 


Materials and methods 
Survey design and methods 


The SEAMAP involves the cooperation of state, federal, 
and university partners for the collection of fishery- 
independent data. The SEAMAP-Gulf of Mexico ground- 
fish trawl survey is conducted twice each year during 
summer and fall at distances from shore up to 370 km (200 
nautical miles) and at depths up to 109.7 m (60 fathoms). 
The summer survey typically occurs in June and July but 
can start as early as May and end as late as August. The 
fall survey typically occurs during October and November 
but can start as early as September and end as late as 
December. Sampling is done by using a stratified random 
design and a 12.2-m (40-ft) trawl net with stretched mesh 
of 40.3-41.9 mm. Environmental conditions, including 
salinity and temperature, were measured either at 3 posi- 
tions in the water column, the surface, mid-water, and bot- 
tom, or through a profile of the full water column at each 
sampling location. 

Prior to fall 2008, sampling was stratified by depth 
(with 23 depth strata, including 1-fathom strata for 5-20 
fathoms, a 2-fathom stratum of 20—22 fathoms, a 3-fathom 
stratum of 22—25 fathoms, 5-fathom strata for 25-50 fath- 
oms, and a single stratum of 25-50 fathoms), by geographic 
area (based on groupings of 2-3 shrimp statistical zones 
[SSZ], which are defined by the National Marine Fisheries 
Service), and by time of day (day or night) (Nichols”). Trawl 
tows were conducted perpendicular to the depth contours 
and covered the entire depth stratum at each sampling 
location. Tow time ranged from 10 to 55 min; if sufficient 
depth coverage was not possible in a single 55-min tow, 
multiple consecutive tows were conducted. Sampling was 
limited to SSZ 11-21. A survey design change occurred 
between summer and fall 2008. Under the new survey 
design, selection of sampling locations was stratified by 
individual SSZ (Fig. 1), and sampling extended into the 
eastern Gulf of Mexico. The stratifications for time of day 
and for depth were dropped. In 2014, depth stratification 
resumed with 2 depth strata: 5-20 fathoms and 20-60 
fathoms (Pollack et al.*). The number of randomly selected 
sampling locations (stations) within each stratum that 
combined geographic area and depth was proportional to 
the geographic area of that stratum. The trawl was towed 
for 30 min at each sampling location. 

For all tows, abundance of each species in the catch was 
recorded. In the case of an extremely large catch (typi- 
cally >22.7 kg), the catch was subsampled and total abun- 
dance of each species was extrapolated. For blue crab, 
carapace widths (CWs), measured as the distance between 
the tips of the lateral spines, of up to 20 individuals 


2 Nichols, S. 2004. Derivation of red snapper time series from 
SEAMAP and groundfish trawl] surveys. Southeast Data Assess. 
Rev. SEDAR7-DW-01, 26 p. [Available from website.] 

3 Pollack, A. G., D. S. Hanisko, and G. W. Ingram Jr. 2016. 
Wenchman abundance indices from SEAMAP groundfish sur- 
veys in the northern Gulf of Mexico. Southeast Data Assess. Rev. 
SEDAR49-DW-19, 27 p. [Available from website.] 
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Figure 1 


A map of the study area in the northwestern Gulf of Mexico showing the locations 
where blue crab (Callinectes sapidus) were collected in groundfish trawl surveys as 
part of the Southeast Area Monitoring and Assessment Program (SEAMAP). Small 
gray dots indicate locations of trawl tows for the study period 2000-2019. Black circles 
indicate locations where adult female blue crab were caught, for reproductive analyses, 
during the SEAMAP groundfish surveys conducted in 2017, with circle size propor- 
tional to the number of blue crab collected from that location. The boxes outlined in 
gray and the numerals inside them indicate shrimp statistical zones (defined by the 
National Marine Fisheries Service), which were used in the stratified random sampling 
design. The black line indicates the 109.7-m (60-fathom) depth contour, which is used 
as the outer limit for SEAMAP groundfish surveys. 


(chosen haphazardly from the total catch) and the sex of 
every fifth measured crab were recorded for each station 
(GSMFC%). 


Long-term patterns of distribution and relative abundance 


Our analyses focused on the period 2000-2019 and 
included only offshore areas that were sampled before 
and after the change in survey design in 2008 (SSZ 11 and 
SSZ 13-21). Chandeleur Sound in Louisiana (SSZ 12) was 
excluded because this inshore area was not consistently 
sampled. Over this 20-year period, 7-119 trawl tows were 
conducted in each SSZ per year (Suppl. Tables 1 and 2) 
(online only). Data collected prior to 2000 were not included 
because of uncertainty surrounding species identifica- 
tion. Data were filtered to include only tows of 12.2-m 
trawl nets with stretched mesh sizes ranging from 40.3 
to 41.9 mm and conducted at depths >9.1 m (>5 fathoms). 
Any sampling events with incomplete data for geographic 
location or vessel speed or with operation codes indicat- 
ing a bad tow (e.g., if the doors came up crossed, if the 


* GSMFC (Gulf States Marine Fisheries Commission). 2016. 
SEAMAP operations manual for trawl and plankton surveys, 61 
p. GSMFC, Ocean Springs, MS. [Available from website.] 


codend came untied, or if the net was fouled or torn) were 
excluded from analysis. 

Maturity stage is not recorded for blue crab; therefore, it 
was necessary to classify crab as adult or juvenile on the 
basis of CW. We classified blue crab >125 mm CW as adults, 
following the cutoff used in recent Gulf-wide (VanderKooy, 
2013) and state-level (West et al.°) blue crab stock assess- 
ments. Because not all blue crab are measured and only 
a small proportion are assigned a sex, it was necessary to 
estimate the catch of adult female blue crab from each trawl 
tow by using the total catch multiplied by the proportion of 
measured blue crab from the tow that were >125 mm CW 
and by the proportion of blue crab from the tow that had 
their sex determined and that were classified as female: 


adult female catch = total catch x proportion adult x 
proportion female. 


For a small number of trawl tows (n=385 tows, 4.3% 
of the tows in the data set used in analysis), sex was 
not recorded, even though adult blue crab were caught 
during the tow. For these tows, the mean proportion of 


5 West, J., H. Blanchet, and P. Cagle. 2019. Update assessment of 
blue crab in Louisiana waters: 2019 report, 32 p. Louisiana Dep. 
Wildl. Fish., Baton Rouge, LA. [Available from website.] 
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crab that were female for that year was used. Sex was 
not determined for blue crab in 2004; therefore, the over- 
all mean proportion that was female (0.930) was used. 
Catch per unit of effort (CPUE) of adult female blue crab 
was calculated by dividing the adult female catch (mea- 
sured as number of individuals) by the area swept by 
the trawl (calculated by using the recorded speed and 
duration of a tow and by assuming a 12.2-m net spread). 
Values of CPUE are given as number of individuals per 
hectare. 

To examine overall temporal and spatial patterns in 
CPUE of blue crab in the SEAMAP survey area, adult 
female catch was modeled by using a 2D generalized 
additive model (Potts and Rose, 2018) with a negative 
binomial error distribution and log link function. Mul- 
tiple candidate models were fit and compared by using 
Akaike information criteria (AIC) to determine the most 
appropriate model. Adult female catch was the response 
variable. Year was included in all models as a categori- 
cal variable. Geographic location was included as a 2D 
smooth function with thin-plate spline bases. Season 
(summer or fall) was included as an additional temporal 
covariate. All models also included an offset term for the 
area swept by the trawl to account for differences in sam- 
pling effort among trawl] tows. Following model selection, 
a Wald chi-square test was used to test for significance of 
predictor variables. 

An index of abundance for the spawning stock in the 
northwestern Gulf of Mexico was estimated to allow 
examination of population trends, accounting for spatial 
and temporal variation in sampling effort. Estimated mar- 
ginal means for each year were extracted from the model 
output, reflecting the main effect of year averaged over the 
full study area and across both seasons. Estimated mar- 
ginal means were then standardized relative to the aver- 
age annual estimated marginal mean across the study 
period by using a Z-transformation, providing an index of 
abundance in units of standard deviations above or below 
the mean for the 20-year period for which we analyzed 
survey data. 


Relative importance of offshore and estuarine waters 


To assess the importance of offshore spawning areas rela- 
tive to estuarine spawning areas for the spawning stock of 
blue crab in the Gulf of Mexico, total abundance of adult 
female blue crab was estimated and compared between off- 
shore and estuarine waters of the study area. Abundance 
estimates were limited to the areas west of the Mississippi 
River (SSZ 13-21), because of the lack of consistent and 
spatially widespread fishery-independent trawl sampling 
for blue crab in estuarine waters of the state of Mississippi. 
Temporal coverage of abundance estimates was limited to 
the period 2006—2018 because of limitations on availabil- 
ity of data for sampling in estuarine areas in other years. 
Observed CPUE was used as a proxy for density of blue 
crab. Because gear efficiencies are unknown for the trawl 
surveys examined in our study, we were unable to cor- 
rect for gear efficiency; actual densities of blue crab, and 


therefore abundances, are likely substantially higher than 
estimates calculated in our study. 

Offshore abundance of blue crab was estimated by using 
the annual mean observed CPUE of adult female blue 
crab from the summer SEAMAP groundfish trawl survey 
and the total area of potential offshore spawning habitat, 
defined as all waters in the Gulf of Mexico west of the Mis- 
sissippi River (SSZ 13-21) from the estuarine mouth or 
barrier islands to the 109.7-m (60-fathom) depth contour 
(depth contours were derived from gridded bathymetry 
data from the General Bathymetric Chart of the Oceans 
(GEBCO_2020 Grid, available from website). Annual mean 
CPUE values, meaning numbers of individuals per hect- 
are, were multiplied by the total area of potential offshore 
spawning habitat to estimate total abundance of blue crab 
in offshore waters for each year. These estimates of off- 
shore abundance were compared with similar estimates of 
abundance in estuarine waters of Texas and Louisiana, by 
using CPUE from fishery-independent trawl surveys con- 
ducted in estuarine areas by the Texas Parks and Wildlife 
Department (TPWD) and Louisiana Department of Wild- 
life and Fisheries (LDWF). Only data collected during 
June and July were included in these analyses because 
these summer months represent the time of peak adult 
female abundance in estuarine waters of Louisiana and 
Texas and also correspond with the timing of the summer 
SEAMAP trawl! survey. 

The TPWD conducts fishery-independent sampling by 
using 6.1-m trawl nets with 38.1-mm stretched mesh, 
towed for 10 min at a speed of 4.8 km/h. Sampling is based 
on a stratified random design, with each of the 8 major 
bay systems in Texas (Sabine, Galveston, Matagorda, San 
Antonio, Aransas, Corpus Christi, Upper Laguna Madre, 
and Lower Laguna Madre) divided into grid cells that 
were 1.85 by 1.85 km. Each month, 10-20 grid cells in 
each major bay system are randomly selected for trawl 
sampling, with the number of cells depending on the bay 
system. Statewide annual mean CPUE of adult (=125 mm 
CW) females was calculated for the period June—July 
and then was multiplied by the total available habitat 
area to estimate the total abundance of mature females 
in estuarine waters of Texas for each year. Habitat area 
was defined as the total area of the Texas coastal bay sys- 
tems extracted from geographic information systems data 
obtained from TPWD (available from website). 

The LDWF conducts fishery-independent sampling 
by using 4.9-m flat trawl nets with a 38.1-mm stretched 
mesh body and a 12.7-mm stretched mesh bag, towed for 
10 min. Sampling occurs at fixed stations, semimonthly 
during November—February and weekly during March— 
October. Vessel speed is not recorded but typically ranges 
from 2.6 to 4.8 km/h (Cagle®). For CPUE calculations, an 
average vessel speed of 3.7 km/h was assumed because it 
is the midpoint of the typical range. Because stations were 
added and deleted during the period of interest, we used 
only data from stations that were sampled in 2018 and had 


® Cagle, P. 2020. Personal commun. Louisiana Dep. Wildl. Fish., 
P.O. Box 98000, Baton Rouge, LA 70898. 
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been sampled for at least 8 years during the period 2006-— 
2018. Additionally, only stations west of the Mississippi 
River were included, to correspond with the spatial extent 
of the SEAMAP data used for offshore abundance esti- 
mates. Annual mean CPUE of adult (=125 mm CW) female 
blue crab was then multiplied by the total available habi- 
tat area to estimate the total abundance of mature females 
in estuarine waters of Louisiana west of the Mississippi 
River for each year. Total available estuarine habitat area 
in Louisiana was calculated by using the NOS80K shore- 
line data set from the NOAA Office of Ocean Resources 
Conservation and Assessment (NOAA’). This shoreline 
layer was created in 1994, however, and substantial land 
loss has occurred in Louisiana since that time (Couvillion 
et al.°). The total available estuarine habitat in the 
Louisiana portion of their study area was, therefore, cor- 
rected on the basis of observed rates of coastal land loss in 
our study area (Couvillion et al.°). 


Condition and reproductive output of spawning females 


Mature female blue crab were collected offshore during 
the summer and fall SEAMAP groundfish trawl surveys 
in 2017 for analyses of overall condition and reproductive 
output. Mature females were haphazardly selected from 
the haul of each trawl tow, individually bagged and frozen, 
and sent to the Gulf Coast Research Laboratory in Ocean 
Springs, Mississippi, for analysis. During the summer sur- 
vey (30 May—2 July 2017), 109 blue crab were collected 
from 48 stations, and during the fall survey (13 October—8 
November 2017), 11 females were collected from 10 sta- 
tions (Fig. 1). 

Following thawing, CW and total mass were recorded for 
each crab. Gills were inspected for the presence of the par- 
asitic gooseneck barnacle Octolasmis muelleri and for the 
parasitic nemertean worm Carcinonemertes carcinophila. 
Spermathecae, the paired sperm storage organs, were 
located to determine if the female had recently mated. A 
female was determined to have recently mated (within 
the past 3-5 weeks; Wolcott et al., 2005) if her spermathe- 
cae had a large size, firm texture, and pink color (Wolcott 
et al., 2005; Ogburn and Habegger, 2015). Because females 
mate only once in their lifetime, following their terminal 
molt, large, firm, pink spermathecae serve as an indica- 
tion that a female moved offshore soon after her mating 
or that a female mated while offshore. Each female was 
determined to be spawning capable if she had developed, 
full ovaries or carried an egg mass. 

For ovigerous females, the egg mass was carefully 
separated from the pleopods and weighed. A subsam- 
ple <1 g (typically 500-1000 eggs) was then removed 


7 NOAA. 1994. NOS80K: medium resolution digital vector U.S. 
shoreline shapefile. Strategic Environ. Assess. Div., Off. Ocean 
Resour. Conserv., Natl. Ocean Serv., NOAA, Silver Spring, MD. 
[Available from website, accessed May 2020.] 

e Couvillion, B. R., H. Beck, D. Schoolmaster, and M. Fischer. 
2017. Land area change in coastal Louisiana 1932 to 2016: 
U.S. Geological Survey scientific investigations map 3381, 16 p. 
[Pamphlet.] [Available from website, accessed May 2020.] 


and weighed, and the total number of eggs in the sub- 
sample was counted by using a stereoscopic microscope. 
Total fecundity was estimated for each ovigerous female 
by dividing the number of eggs in the subsample by the 
mass of the subsample and multiplying by the mass of 
the entire egg mass. Relationships between fecundity 
and body size were examined by using generalized lin- 
ear models (GLMs) with a Gaussian distribution and log 
link function and with fecundity as the response variable. 
Separate GLMs were fit to model fecundity as a function 
of CW and body mass (excluding the mass of the egg 
mass). Because of missing data for body mass of 16 of 
the 73 ovigerous blue crab, it was not possible to directly 
compare the model fits by using AIC or the AIC corrected 
for bias of small sample sizes (AICc) because the 2 mod- 
els used slightly different data sets. For comparisons, the 
GLM with CW as the predictor was refit by using only 
individuals with data for both CW and body mass, and 
then the models were compared by using AICc. 

The embryonic developmental stage for each egg sam- 
ple was determined by visual examination under a com- 
pound microscope and classified as early (stages 1-3 of 
DeVries et al., 1983), middle (stages 4-6 of DeVries et al., 
1983), and late (stages 7—9 of DeVries et al., 1983). The 
percentage of embryos that were developing normally 
(Darnell et al., 2009) was visually assessed for each egg 
mass on the basis of 20 eggs chosen at random from the 
subsample. Egg size for 20 normally developing embryos 
was determined by measuring along the long axis and 
the perpendicular axis with a compound microscope. Egg 
volume was calculated for each measured egg as the vol- 
ume of a prolate spheroid by using the short and long 
axis measurements. Egg volume measurements were 
averaged for each egg mass, and this mean volume for 
each egg mass was then used in subsequent analyses. 
Egg volume was analyzed by using a GLM with a gauss- 
ian distribution and a log link function and with egg vol- 
ume as the response variable. Egg volume was modeled 
as a function of CW and embyronic developmental stage 
(early, middle, and late). 


Results 
Long-term patterns of distribution and relative abundance 


A total of 8996 trawl tows were conducted in offshore 
waters across the study area in 2000-2019. During that 
time, 97.7% of summer sampling occurred in June and 
July, with 1.3% occurring in May and 1.0% occurring in 
August. Fall sampling was primarily (97.7%) conducted 
during October and November, with 2.2% occurring in 
September and 0.1% occurring in December. Depths sam- 
pled were 9.2—109.5 m. Salinities measured at the bottom 
of the water column at sampling locations were 10.9-39.5, 
with 97.3% of trawl tows conducted in waters with salini- 
ties >30.0. Bottom temperatures were 16.0—30.5°C. 

Blue crab were caught in 19.4% of all trawl tows, with 
11,021 caught during 2000-2019 across the full study 
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area. Recorded CWs were 3—281 mm (median: 136 mm), 
with 99.5% of blue crab having CWs of 20-200 mm. 
Although sex was not recorded for all blue crab, 93.0% 
of the 1597 blue crab for which sex was determined were 
female. Adult female blue crab were caught in 17.1% of all 
trawl tows, with a total of 6355 caught. The overall mean 
CPUE of adult females was 0.25 individuals/ha (standard 
error of the mean [SE] 0.01) over the full study period, with 
annual mean CPUE ranging from 0.06 individuals/ha 


(SE 0.02) in 2013 to 0.53 individuals/ha (SE 0.09) in 2010.: 


Annual frequencies of occurrence for adult female blue 
crab ranged from 8.5% in 2000 to 25.7% in 2011 (Table 1). 
The maximum annual mean CPUE of adult females 
observed during the summer surveys was 0.87 individu- 
als/ha (SE 0.31) in 2006, and the maximum annual mean 
CPUE of adult females observed during the fall surveys 
was 0.29 individuals/ha (SE 0.09) in 2011. 

Following comparisons of candidate models, the selected 
model included year, season, an interaction of year and 
season, and separate smooth terms for geographic location 
in each season (Table 2), with the selected model explain- 
ing 47.9% of the deviance in catch. Results of the Wald 


Table 2 


Candidate models used to estimate catch of adult female 
blue crab (Callinectes sapidus) collected from the north- 
west Gulf of Mexico between 2000 and 2019. Model selec- 
tion criteria also are given for each model: the difference in 
Akaike information criterion (AIC) between each model and 
the model with the lowest AIC (AAIC), AIC weight (wajc), 
and the proportion of deviance explained. The characters 
“s()” indicate that a continuous variable was included as 
a smooth term. All models included an offset term for the 
area swept by the trawl to account for differences in sam- 
pling effort among trawl tows. 


Model formulation AIC AAIC Waic D? 


Year + season + 12,947.7 0.00 1 
(year x season) + 
s(geographic 
location) x season 
Year + season + 13,072.6 
s(geographic 
location) x season 
Year + season + 13,302.3 
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Table 1 


Summary information by year for adult female blue crab 
(Callinectes sapidus) collected in the northwestern Gulf 
of Mexico during groundfish trawl surveys conducted in 
2000-2019 as part of the Southeast Area Monitoring and 
Assessment Program. The number of trawl tows; total 
catch (number of individuals); mean catch per unit of 
effort (CPUE), calculated by dividing catch by the area 
swept by the trawl (number of individuals per hectare) 
and given with standard error of the mean (SE) in paren- 
theses; and percent occurrence of female adult blue crab 
are provided. 


No. of Total Mean Percent 
Year trawl tows catch CPUE (SE) occurrence 


2000 585 159 
2001 525 179 
2002 575 
2003 528 
2004 533 
2005 517 
2006 529 
2007 510 
2008 642 
2009 689 
2010 431 
2011 381 
2012 336 
2013 261 
2014 342 
2015 346 
2016 309 


0.14 (0.06) 8.55 
0.14 (0.03) 10.86 
0.19 (0.03) 17.39 
0.26 (0.05) 16.48 
0.26 (0.05) 14.63 
0.31 (0.04) 22.63 
0.48 (0.16) 18.34 
0.26 (0.05) 16.86 
0.23 (0.05) 13.86 
0.33 (0.05) 20.90 
0.53 (0.09) 23.43 
0.49 (0.08) 25.72 
0.16 (0.04) 15.18 
0.06 (0.02) 9.58 
0.11 (0.02) 15.79 
0.12 (0.02) 18.21 
0.16 (0.03) 20.71 
2017 333 0.26 (0.05) 19.82 
2018 308 0.13 (0.02) 17.86 
2019 0.17 (0.03) 18.35 


(year x season) + 
s(geographic 
location) 
Year + season + 13,410.7 
s(geographic 
location) 
Year + s(geographic 13,636.8 689.07 
location) 
Year + season 14,560.2 1612.55 
Year 14,797.0 1849.35 


chi-square test indicate that all predictors had a significant 
effect on catch (Table 3). Averaged over the full study period, 
CPUE of adult female blue crab was higher during the sum- 
mer surveys (0.39 individuals/ha [SE 0.03]) than during the 
fall surveys (0.11 individuals/ha [SE 0.01]), although this 
pattern was not consistent across all years (Fig. 2). The 
spatial distribution of adult female blue crab also differed 
between the 2 seasons (Table 3). In summer surveys, CPUE 
was generally highest in nearshore regions off Louisiana 
and along the central and upper coasts of Texas, primarily 
at depths <36.6 m (<20 fathoms). The CPUE from fall sur- 
veys was generally highest off southeastern Louisiana and 
at intermediate distances from shore (Fig. 3). 

Estimated marginal mean CPUE for each year, stan- 
dardized relative to the 20-year average, provides a stan- 
dardized index of abundance for the offshore spawning 
stock of blue crab in the northwestern Gulf of Mexico 
(Fig. 4). With the exception of 2010 and 2011, the index 
of abundance for each year was within 1.5 standard devi- 
ations of the mean value. The index indicates a trend of 
slight increases in abundance from 2000 to 2009, followed 
by a large increase in 2010, a sustained high level of abun- 
dance in 2011, and a marked decrease after 2011. After 
a low in 2012, the index has steadily increased back to 
average levels of abundance. 
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Table 3 


Results from the Wald chi-square (x) test for significance 
of predictor variables in the 2D generalized additive model 
used to examine overall temporal and spatial patterns in 
catch of adult female blue crab (Callinectes sapidus) col- 
lected during trawl surveys in the northwestern Gulf of 
Mexico between 2000 and 2019. The characters “s()” indi- 
cate that a continuous variable was included as a smooth 
term. Degrees of freedom (df) and estimated degrees of 
freedom (edf) are shown for parametric factor terms and 
smooth terms, respectively. 


Source df or edf "2 P 


Year 19 240.40 <0.0001 

Season 1 21.54 <0.0001 

Year x season 19 174.71 <0.0001 

s(geographic location): 66.16 745.00 <0.0001 
summer 

s(geographic location): fall 69.20 530.10 <0.0001 


Relative importance of offshore and estuarine waters 


During 2006-2018, considering only SSZ 13-21, CPUE 
of adult females from summer surveys in offshore waters 
of the Gulf of Mexico averaged 0.38 individuals/ha 
(SE 0.03), compared with the CPUE values of 0.46 
individuals/ha (SE 0.03) in estuaries of Texas and of 3.92 
individuals/ha (SE 0.17) in estuaries of 
Louisiana. Total available habitat area 
was estimated to be 11,473,937.0 ha 
for the offshore SSZ 13-21 in the 
SEAMAP study area, 578,310.7 ha for 
estuaries of Texas, and 853,436.2 ha 
for estuaries of Louisiana west of the 
Mississippi River. Rough calculations 
made by using CPUE of adult females 
(as a proxy for density) and total habi- 
tat area indicate that estuarine areas 
during June—July supported a min- 
imum abundance of 1,248,675 adult 
female blue crab in 2017 and a maxi- 
mum abundance of 8,631,449 adult 
female blue crab in 2006 and that the 
offshore SEAMAP survey area during 
June—July supported a minimum abun- 
dance of 664,203 adult female blue crab 
in 2013 and a maximum abundance 
of 9,509,020 adult female blue crab in 
2010 (Fig. 5). On average, across the 
13-year period of 2006-2018, the off- 
shore SEAMAP survey area supported 
48.8% of the total spawning stock in the 
northwestern Gulf of Mexico, although 
this percentage ranged from 26.5% in 
2014 to 79.7% in 2017. 


CPUE (individuals/ha) 


Condition and reproductive output of spawning females 


Collection of blue crab for reproductive analyses occurred 
on each state and federal leg of the SEAMAP groundfish 
trawl surveys during summer and fall 2017. In total, 
120 mature female blue crab were collected, with CWs 
ranging from 110 to 209 mm. The vast majority of blue 
crab (109 of 120) were collected during the summer 
survey (30 May—2 July 2017), and the rest of the indi- 
viduals (11 of 120) were caught during the fall survey 
(13 October—8 November 2017). 

Of the 109 mature females collected during summer, 
95 individuals (87.2%) were spawning capable, as evidenced 
by an external egg mass (65 individuals) or mature ova- 
ries (30 individuals). All embryonic developmental stages 
were observed in ovigerous blue crab collected during the 
summer survey, with 56.9% (37 individuals) in the early 
developmental stage, 10.8% (7 individuals) in the middle 
developmental stage, and 32.3% (21 individuals) in the late 
developmental stage. Additionally, 26.2% (17 individuals) 
of the ovigerous females collected during the summer sur- 
vey also possessed developed, full ovaries, indicating the 
capacity to produce at least one additional clutch of eggs 
following the clutch carried at the time of collection. Only 
1 female (0.9%) collected during the summer survey had 
mated recently. 

Of the 11 mature females collected during the fall sur- 
vey, 6 individuals (54.4%) were spawning capable, as evi- 
denced by an external egg mass (5 individuals) or 
developed, full ovaries (1 individual). Of these 5 oviger- 
ous individuals, 4 crab carried eggs in the early 


-® Annual 
~- Summer 
-7- Fall 


Figure 2 


Observed mean annual (solid line) and seasonal (broken lines) catch per unit 
of effort (CPUE) for adult female blue crab (Callinectes sapidus) collected 
during groundfish surveys conducted in the northwestern Gulf of Mexico from 
2000 through 2019 as part of the Southeast Area Monitoring and Assessment 
Program. Values of CPUE were calculated by dividing catch by the area swept 
by the trawl. Error bars indicate standard errors of the mean. 
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Figure 3 


Spatial patterns of catch per unit of effort (CPUE) for adult female blue crab (Callinectes sapidus) collected during ground- 
fish trawl surveys conducted in the northwestern Gulf of Mexico in summer (left panels) and fall (right panels) across the 
study period 2000-2019. The top panels show observed CPUE (number of individuals per hectare), binned into 0.25° grid 
cells and averaged across all years on a log scale. White areas represent areas that were not sampled, and gray grid cells 
represent areas with no catch of adult female blue crab over the full study period. The middle panels show interannual 
variability in catch rates, calculated as the coefficient of variation of observed annual mean CPUE. The bottom panels show 
predicted CPUE of adult female blue crab for each season from the generalized additive model, averaged across all years and 
standardized by using a Z-transformation based on the mean and standard deviation for that particular season. Values are 
depicted in units of standard deviations above (positive values) or below (negative values) the mean. White contours delineate 
a standardized value of zero, representing the mean for that season. Black lines represent the 36.6-, 73.2-, and 109.7-m (20-, 
40-, and 60-fathom) depth contours. 


developmental stage and 1 crab carried eggs in the late fecundity. This relationship was evident when both CW 
stage. None of the females collected during the fall sur- (P<0.0001, deviance explained [D?]=0.539) and body mass 
vey had mated recently. (P<0.0001, D?=0.712) were considered (Fig. 6). Results of 

Fecundity was 1.07—8.93 million eggs (mean: 3.27 mil- the comparison of model fits based on AICc indicate that 


lion eggs [SE 0.15 million]) and was strongly related to body mass was a better predictor of fecundity than CW, 
size of female crab, with larger females having higher with an AlCc weight >0.98. 
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Figure 4 


Standardized index of abundance for adult female blue crab (Call- 
inectes sapidus) in the northwestern Gulf of Mexico during the study 
period 2000-2019. The gray shaded area indicates the 95% confi- 
dence interval. The horizontal dotted line indicates a value of zero, 
representing the average annual predicted catch per unit of effort 


across the full study period. 
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Figure 5 


Estimated total abundance (in millions) of adult female blue crab 
(Callinectes sapidus) in offshore (solid line) and estuarine (dashed 
line) areas of the northwestern Gulf of Mexico west of the Missis- 


sippi River during June and July in 2006-2018. 


On average across both seasons, 91.7% (SE 0.01) of eggs 
were fertilized and embryos were developing normally, 
with 91.6% (SE 0.01) of embryos developing normally in 
the summer and 93.0% (SE 0.05) of embryos developing 
normally in the fall. Eggs that were not developing nor- 
mally were either unfertilized or fertilized but had ceased 
development in an earlier developmental stage than the 
majority of the egg mass. Average egg size was 243.6 pm 
(SE 0.5) along the long axis and 236.1 pm (SE 0.5) along 


the perpendicular axis. Average egg volume was 
0.0073 mm? (SE 0.0001). There was a significant 
relationship between egg volume and embryonic 
developmental stage (F=29.76, P<0.0001), with a 
general increasing trend as the embryos develop 
further (Fig. 7). There was no significant relation- 
ship between egg volume and maternal CW 
(F=2.371, P=0.1286) and no significant interac- 
tion between maternal CW and embryonic devel- 
opmental stage (F=2.363, P=0.1023). 

The parasitic barnacle O. muelleri was com- 
monly found on the gills of blue crab, with prev- 
alence of 66.1% in summer and prevalence of 
100.0% in fall. The nemertean worm C. carci- 
nophila that is a predator of crab eggs was less 
common, with prevalence of 6.4% in summer and 
prevalence of 0.0% in fall. 


Discussion 


Results of this study indicate that the offshore 
waters of the northwestern Gulf of Mexico are 
an important spawning habitat for blue crab. 
Mature female blue crab are routinely caught 
in offshore waters of the northwestern Gulf of 
Mexico each year in the SEAMAP groundfish 
trawl survey, as far as 200 km from shore (Fig. 3). 
Our results support the substantial numbers of 
spawning blue crab in offshore areas of the Gulf 
of Mexico documented in previous studies (Gelpi 
et al., 2009, 2013). 

The CPUE of adult female blue crab in offshore 
waters was greater in the summer than in the 
fall during most years (Fig. 3). This difference in 
CPUE between summer and fall may be indica- 
tive of a postspawning return to inshore habitats, 
as proposed by Tankersley et al. (1998). Given 
that findings from more recent studies indicate 
that migration back to estuaries is unlikely 
(Hench et al., 2004; Darnell et al., 2012), this 
pattern more likely is indicative of an offshore 
population decline from natural mortality due 
to predation and the high energy expenditures 
of migration and reproduction. It is important to 
note that the SEAMAP-Gulf of Mexico groundfish 
surveys are limited to summer (primarily June— 
July) and fall (primarily October-November). 
Blue crab recruitment in the northwestern Gulf 
of Mexico generally peaks from August to October 
(Rabalais et al., 1995; Perry et al. 2003), which corresponds 
with the observed peak in adult female densities during 
the summer SEAMAP trawl survey. Given that earlier 
recruitment peaks are observed in some areas (Rabalais 
et al., 1995; Sullivan and Neigel, 2017), late spring may 
also be an important time for offshore spawning in this 
region (data are not available to test this hypothesis). 

Throughout the study period, the highest densities of 
adult female blue crab were present in the waters off 
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Figure 6 

Relationships between fecundity and (A) carapace width and (B) 
body mass for ovigerous female blue crab (Callinectes sapidus) col- 
lected during the groundfish trawl survey conducted in 2017 in the 
northwestern Gulf of Mexico as part of the Southeast Area Monitor- 
ing and Assessment Program. Lines indicate relationships from gen- 
eralized linear models. Shaded regions represent the 95% confidence 
interval. D’=deviance explained. 


Louisiana, and this pattern was consistent for both sum- 
mer and fall surveys (Fig. 4). This spatial pattern reflects 
reported inshore fishery landings, with areas off Louisi- 
ana accounting for 88% of the blue crab landings in the 
Gulf of Mexico in 2018 (NMFS’). The large numbers of 
blue crab caught offshore of Louisiana relative to catch 
in other regions likely reflect the low salinities present in 
estuaries of Louisiana. Egg hatching and larval survivor- 
ship is negatively affected by salinities <20.1 (Costlow 
and Bookhout, 1959); during the spawning migration, 
therefore, females move to spawning grounds suitable 
for offspring survival. It is possible that the low salini- 
ties in estuaries of Louisiana result in greater numbers 


° NMFS (National Marine Fisheries Service). 2018. Commercial 
landings query, release 3.1.0.2. Off. Sci. Technol., Natl. Mar. 
Fish. Serv, Silver Spring, MD. [Available from website, accessed 
March 2020.] 


of females migrating out of the estuaries and 
continuing to migrate into the offshore environ- 
ment, compared with the numbers of females 
that migrate out of estuaries in other regions. 
Although CPUE was highest off Louisiana in 
both seasons, the spatial distribution of the 
CPUE differed between summer and fall (Fig. 4). 
During summer, adult females were primarily 
caught in nearshore regions, at depths <36.6 m 
(<20 fathoms). During fall, adult females were 
primarily caught at depths of 36.6—73.2 m (20-40 
fathoms) off both Texas and Louisiana. 

Ogburn and Habegger (2015) conducted a sim- 
ilar examination of blue crab spawning in off- 
shore waters of North Carolina, South Carolina, 
Georgia, and Florida, by using data from the 
SEAMAP-South Atlantic trawl surveys, with sim- 
ilar results. The highest CPUE of mature female 
blue crab occurred during the summer, females 
were actively spawning offshore, and a substantial 
portion of the spawning stock in Georgia occurred 
in offshore waters. Using only data from summer 
surveys, we observed a maximum annual mean 
CPUE of 0.87 individuals/ha (SE 0.31) (in 2006), 
a level substantially lower than the maximum of 
1.58 individuals/ha observed (in 1991) by Ogburn 
and Habegger (2015). When interpreting this 
difference in CPUE, however, it is important to 
consider the differences in survey design between 
the SEAMAP trawl surveys conducted along the 
Atlantic coast of the southeastern United States 
and in the Gulf of Mexico. The SEAMAP-South 
Atlantic survey samples at depths of 4.0-10.0 m 
(2-5 fathoms), and the SEAMAP-Gulf of Mexico 
survey samples at depths of 9.1-109.7 m (5-60 
fathoms). If the SEAMAP-Gulf of Mexico data set 
examined in this study is limited to depths <15 m, 
mean adult female CPUE was 1.96 individuals/ 
ha (SE 0.97) in summer 2006, 1.58 individuals/ha 
(SE 1.13) in summer 2010, and 1.53 individuals/ha 
(SE 1.09) in summer 2000. 

Although densities of adult female blue crab observed 
in the SEAMAP study area were relatively low, the off- 
shore areas sampled in our study support a substantial 
portion (up to 79.7%, depending on the year) of the 
spawning stock in the northwestern Gulf of Mexico. We 
acknowledge that our estimates of total abundance of the 
spawning stock are rough, with multiple potential 
sources of error. A size-based cutoff was used to classify 
blue crab as adults or juveniles, rather than a visual mor- 
phological examination. Although the fishery-indepen- 
dent trawl surveys conducted by the LDWF and TPWD 
record maturity stage for blue crab, this information is 
not recorded in the SEAMAP trawl surveys. For consis- 
tency, we therefore classified all blue crab >125 mm CW 
as adults for data from all 3 survey programs. Addition- 
ally, these estimates do not account for differences in 
gear or gear efficiencies in the 3 survey programs. Actual 
abundances, therefore, are likely to be substantially 
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Figure 7 


Mean egg volume for each embryonic developmental stage 
for ovigerous female blue crab (Callinectes sapidus) col- 
lected during the groundfish trawl survey conducted in 
2017 in the northwestern Gulf of Mexico as part of the 
Southeast Area Monitoring and Assessment Program. 
Mean egg volumes are estimated marginal means from 
the generalized linear model. Error bars indicate the 
back-transformed 95% confidence intervals. The different 
letters above the bars indicate statistically significant dif- 
ferences between developmental stages. 


higher than the estimates presented here. It is also 
important to consider the fate of larvae spawned offshore. 
There may be some distance offshore beyond which it is 
highly unlikely that a blue crab larva could successfully 
make it back to an estuary. To date, most larval transport 
modeling efforts for blue crab in the Gulf of Mexico have 
focused on spawning occurring at or near the estuarine 
mouth (e.g., Johnson and Perry, 1999; Jones et al., 2015; 
Criales et al., 2019). Further work is needed to better 
understand the fate of larval blue crab that have been 
spawned offshore. 

Our analyses of the reproductive health of female blue 
crab collected in offshore habitats gave no indication 
that these females are of poor health or that they are 
of any lesser value to the spawning stock than females 
spawning in estuarine waters. Estimated fecundity was 
1.07-8.93 million eggs, with an average of 3.27 million 
eggs (SE 0.15 million). In comparison, Graham et al. 
(2012) reported fecundity estimates of ~0.9—7.0 million 
eggs for ovigerous blue crab collected from estuarine 
waters of Mississippi Sound, and Prager et al. (1990) 
reported fecundity estimates of ~0.5-8.0 million eggs 
(mean: 3.2 million eggs) for ovigerous blue crab col- 
lected from Chesapeake Bay. Egg size in our study was 
also similar to that found in previous studies of oviger- 
ous females collected in estuarine waters. We observed 
egg measurements of 243.6 pm (SE 0.5) along the long 
axis and of 236.1 pm (SE 0.5) along the perpendicular 


axis. In comparison, Graham et al. (2012) reported egg 
diameters of ~215-310 pm (~215-275 pm during the 
summer and fall only) for blue crab from Mississippi 
Sound, and Darnell et al. (2009) reported an average 
egg diameter of 267.5 pm (SE 1.9) for blue crab from 
North Carolina. 

Fecundity of blue crab declines with age, with the 
greatest fecundity in the first clutch of eggs produced by 
a female (Darnell et al., 2009). That females collected in 
offshore and estuarine waters produce similar quantities 
of eggs indicates that offshore females are at a similar 
stage of their reproductive period as females spawning in 
estuarine areas. This notion is supported by our observa- 
tion of a high percentage of eggs that were fertilized with 
embryos developing normally (91.7% [SE 0.01]). Darnell 
et al. (2009) observed a decrease in the percentage of 
embryos developing normally in later clutches, from 
96.7% (SE 0.9) in a female’s first clutch to 55.0% (SE 20.1) 
in its fourth clutch. Assuming a similar pattern in blue 
crab from the Gulf of Mexico, this finding indicates that 
the females assessed in our study were producing their 
first or second clutch of eggs. We observed full, developed 
ovaries in >50% of females collected offshore during the 
summer, indicating that these blue crab would spawn at 
least once more in their lifetime. Our results also indicate 
that oviposition and continued spawning is occurring in 
offshore waters, and we observed all stages of embryonic 
development. 

The results of this study indicate that the coastal and 
offshore waters of the northwestern Gulf of Mexico are 
an important spawning habitat for blue crab. Although 
densities of blue crab offshore were low, these crab were 
primarily adult females and were actively spawning. 
Given the large amount of available shelf habitat, our 
results indicate that offshore areas support a substan- 
tial portion of the spawning stock in the northwestern 
Gulf of Mexico. Fecundity of females caught offshore was 
similar to that reported for estuarine spawning crab, and 
these females are continuing to spawn multiple clutches 
of eggs in offshore waters. Additional research is needed 
to understand the causes of temporal and spatial varia- 
tion in densities of spawning female blue crab in the Gulf 
of Mexico. Furthermore, it is critical that this previously 
understudied portion of the spawning stock of blue crab 
in the Gulf of Mexico be considered and accounted for in 
future management and assessment efforts in the Gulf 
of Mexico. 
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Abstract—Median length at maturity 
was determined for Pacific angel sharks 
(Squatina californica) captured inciden- 
tally in the artisanal elasmobranch fish- 
ery in the Gulf of California in Mexico. 
A total of 306 Pacific angel sharks (192 
females and 114 males) were analyzed. 
The size of the females ranged between 
23 and 100 cm total length (TL), 
whereas males ranged between 25 and 
99 cm TL. The maturity stages of both 
females and males were determined by 
using the development of internal and 
external organs. The results of analysis 
of covariance reveal a significant effect 
of sex (P<0.001) on length at maturity. A 
binary logistic regression was applied for 
each sex to estimate the length at which 
half of the individuals were consid- 
ered mature (L;)). For females, L;) was 
74.41 cm TL (95% confidence interval 
[CI]: 72.81—76.00 cm TL), and the steep- 
ness coefficient (©) was 1.43 (95% CI: 
0.84-2.42). For males, L;y was 77.82 cm 
TL (95% CI: 75.66-79.97 cm TL), and 
the ® was 3.08 (95% CI: 1.98-4.77). 
The maturity ogive produced by using a 
logistic function is suitable to represent 
the development pattern of this species 
because it describes maturity as a grad- 
ual process. 
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The Pacific angel shark (Squatina 
californica) belongs to the family Squa- 
tinidae and inhabits the eastern Pacific 
Ocean. Its geographical distribution is 
discontinuous: from southern Alaska 
to the Gulf of California (GOC) in Mexico 
and from Ecuador to Chile (Compagno 
et al., 2005). This species is a noctur- 
nal opportunistic predator of demersal 
fish species (Pittenger, 1984; Fouts and 
Nelson, 1999). In California, this spe- 
cies has been found in shallow waters 
at depths of 3-46 m, whereas in the 
GOC it has been observed at a depth 
of 183 m. This species has been found 
buried in muddy and sandy bottoms. 
It also has been observed around rocks 
and in kelp forests (Compagno, 1984). 
The demography of the northern pop- 
ulation off California has been stud- 
ied, and the low productivity of this 
population has been modeled (Cailliet 
et al., 1992). The low productivity of 
this species is consistent with the char- 
acteristics of elasmobranchs that make 
them susceptible to overfishing: late 
maturity, low fecundity, and long life 
(Walker, 1998). 

The capture of Pacific angel sharks 
is mostly incidental in the artisanal 
elasmobranch fishery and the finfish 
trawl fishery both in the GOC and on 
the Pacific coast of the Baja California 


Peninsula (BCP) in Mexico. Biological 
information about the Pacific angel 
shark in Mexico is fragmented and lim- 
ited to a few studies on trophic ecology 
and reproduction (Escobar-Sanchez 
et al., 2006, 2011; Romero-Caicedo 
et al., 2016). Using mitochondrial 
DNA, Ramirez-Amaro et al. (2017) 
reported genetic differences between 
specimens from the GOC and Pacific 
coast of the BCP, validating previous 
observations of life history differences 
between specimens from these areas 
(Villavicencio-Garayzar, 1996). Sim- 
ilarly, Marquez-Farias (2007) found 
maximum lengths and lengths at 
maturity for the population of shov- 
elnose guitarfish (Pseudobatus pro- 
ductus) in the GOC that are different 
from those reported for the popula- 
tion along the Pacific coast of the BCP 
(Villavicencio-Garayzar, 1995). 
Because of its historical exploitation 
rates and vulnerability, the Pacific angel 
shark is classified as near threatened 
on the IUCN Red List of Threatened 
Species (Cailliet et al., 2016). The sig- 
nificant fishing effort in the elasmo- 
branch fishery of the GOC (Bizzarro 
et al., 2007) makes understanding the 
life history and population dynamics of 
this species essential to identification of 
strategies for its future management. 
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The objective of this study was to estimate median length 
at maturity (L;,) of the Pacific angel shark in the eastern 
GOC. Differences in life history characteristics have import- 
ant management implications if they result in different 
population responses to fishing pressure. 


Materials and methods 


The artisanal elasmobranch fisheries of Sonora and Sinaloa, 
Mexico (Fig. 1), operate open-hulled fiberglass boats 
that are 5.5-7.6 m long and have outboard motors with 
55-115 hp. Pacific angel sharks were collected at fishery 
landings during 1998-2005 (for details, see Bizzarro et al., 
2009). Length, sex, and maturity stage of specimens were 
recorded. Total length (TL) was measured in centimeters 
from the tip of the rostrum to the end of the stretched cau- 
dal fin. Length-frequency distributions were constructed to 
examine the structure of the sample. Normality of length 
distributions was verified by using the Kolmogorov— 
Smirnov (K-S) test. If normality data were rejected, the 
nonparametric Mann—Whitney test was used to compare 
mean lengths of each sex. A chi-square test was used to 
compare the observed sex ratio of the entire sample and the 
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Figure 1 


A map of the study area, showing the fishing grounds 
(black areas) in the Gulf of California off Sonora and 
Sinaloa, Mexico, where Pacific angel sharks (Squatina 
californica) were captured between 1998 and 2005. 


ratio of maturity stages (immature to mature) against the 
expected value of 1:1. Reproductive status was assessed 
through macroscopic examination following Clark and von 
Schmidt (1965) and Castro (2000). Mature males had hard- 
ened claspers that could rotate from the base by directing 
the claspers anteriorly on the inner side of the shark. The 
inner clasper length was measured from the point of inser- 
tion of the clasper to the clasper tip. The clasper length 
was plotted against TL to show how the clasper grew as an 
individual approached maturity. Females were considered 
mature when they had ripe oocytes (>20 mm in diameter) 
and a thickened uterus or when they were pregnant. His- 
tograms and box plots of lengths by maturity stage were 
constructed to identify outliers. 

For statistical treatment, the classification of maturity 
stages used 0 for immature individuals and 1 for mature 
individuals. Binary logistic regression between TL and 
maturity data was applied to estimate the regression 
coefficients a and b. The significance of the difference in 
maturity stage between sexes was tested with an analysis 
of covariance by applying a general linear model with a 
significance level of 0.01. With the coefficients a and b, the 
proportion of mature sharks by length was determined by 
using the following logistic function (maturity ogive): 


1 


Pm», = 5 
a Gy Fl 
© 


1+ exe 


where Pm = the proportion of mature fish given the length 
(TL); 
Ls = the median length of the sample when half 
of the individuals had reached maturity 
(L;9=a/b); and 
® = the steepness of the model (®=L,,)/a). 


All the parameters and their 95% confidence intervals 
(CIs) were estimated with Minitab!, vers. 19.0 (Minitab, 
State College, PA). 


Results 


A total of 306 Pacific angel sharks (192 females and 114 
males) were measured (y7=19.9, P<0.05). The mean length 
for females was 72.91 cm TL (standard deviation [SD] 
16.88) with a range of 23-100 cm TL; the mean length for 
males was 80.32 cm TL (SD 13.32) with a range of 
25-99 cm TL (Fig. 2A). For both sexes, the length distribu- 
tions had maximum frequencies in the length class of 
85-90 cm TL, a class that mostly corresponds with mature 
individuals (Fig. 2, B and C). Results of analysis of the 
normality of length distributions with the K-S test indi- 
cate non-normal distributions of length for males (K—-S 
test: D=0.209, P<0.01) and females (K-S test: D=0.174, 
P<0.01). The nonparametric test revealed significant 


' Mention of trade names or commercial companies is for identi- 
fication purposes only and does not imply endorsement by the 
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Figure 2 


Length-frequency distributions (A) by sex and by matu- 
rity stage for (B) female and (C) male Pacific angel sharks 
(Squatina californica) captured from the eastern Gulf of 
California in Mexico between 1998 and 2005. 


differences in mean length by sex for the Pacific angel 
sharks captured in the GOC (Mann—Whitney test: 
U=26,818, P=0.001). The mean length for mature females 
was 84.83 cm TL (SD 4.83) with a range of 75-100 cm TL 
(number of specimens [n]=114). The diameters of the larg- 
est oocytes were 20-60 mm and, in pregnant females, 
ranged between 26 and 50 mm (Fig. 3A). The mean length 
of mature males was 86.18 cm TL (SD 5.54) with a range 
of 73-99 cm TL (n=78); clasper lengths ranged from 10 to 
30 mm (Fig. 3B). Length ranges and maturity condition by 
sex are presented in Table 1 and Figure 4. 

The analysis of covariance from the binary logistic 
regression revealed a significant effect of sex (y7=14.9, 
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Figure 3 

The relationship between total length and (A) maximum 
oocyte diameter and (B) clasper length for female and 
male Pacific angel sharks (Squatina californica), respec- 
tively, captured from the eastern Gulf of California in 
Mexico during 1998-2005. Open circles indicate immature 
individuals, black circles represent mature individuals, 
and open triangles indicate pregnant females. 


df=1, P<0.001). The model explained 75.0% of the vari- 
ance. The contribution of TL to variance was 71.5%, and 
sex explained 3.7% of the variance. The coefficients of the 
binary logistic regression for females (a=—51.95, b=0.698) 
were used in the logistic function to describe the estimated 
proportion of mature fish by length. For females, the L;o 
was 74.41 cm TL (95% CI: 72.81—76.00 cm TL), and the ® 
was 1.43 (95% CI: 0.84-2.42) (Fig. 5A). For males, the 
regression coefficients (a=—25.27, b=0.3247) yielded a L; 
of 77.82 cm TL (95% CI: 75.66-79.97 cm TL) and a ® of 
3.08 (95% CI: 1.98—4.77) (Fig. 5B). 


Discussion 


The results of this study, conducted on the coasts of Sonora 
and Sinaloa, indicate differences between the lengths of 
Pacific angel sharks in the GOC and those of individuals 
along the Pacific coast of the BCP. In general, individuals 
from the Pacific coast of the BCP grow larger than those 
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Table 1 


Length-frequency distributions, by sex and maturity stage, for Pacific angel sharks (Squatina californica) captured in 
the eastern Gulf of California in Mexico between 1998 and 2005. 


Female 


Total 
length (cm) 


Immature Mature Total 
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Figure 4 


Box plot of the relationship between total length and matu- 
rity stage for Pacific angel sharks (Squatina californica) 
captured from the eastern Gulf of California in Mexico 
during 1998-2005. The black dots indicate observed lengths 
by sex and maturity stage. The upper and lower parts of 
each box represent the first and third quartiles (the 25th 
and 75th percentiles), and the horizontal line is the median. 
The vertical lines that extend above and below each box 
indicate minimum and maximum values. 


from the GOC. In addition, a latitudinal gradient in the 
maximum reported sizes of Pacific angel sharks has been 
observed, with sharks being larger on the coast of California 
and decreasing in length toward the southern part of the 
BCP on the Pacific coast. For instance, Miller and Lea (1972) 
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reported a maximum length of 152.4 cm TL for individuals 
from the coast of California, and Natanson and Cailliet 
(1986) reported an L;, of 107 cm TL for the same area (~70% 
of the maximum length). For individuals on the northern 
part of the BCP, Cartamil et al. (2011) reported lengths of 
120-130 cm TL; whereas, Villavicencio-Garayzar (1996) 
reported a maximum length of 117 cm TL for individuals on 
the southern part of the BCP. Ramirez-Amaro et al. (2017) 
reported an unusually large size of 180 cm TL for individu- 
als in this same area. For Pacific angel sharks in the Bay of 
La Paz, part of the GOC, Romero-Caicedo et al. (2016) 
reported that the largest sizes were 90-100 cm TL and esti- 
mated that the L;, for males and females was 75.6 and 
77.7 cm TL, respectively (~77% of the maximum length). In 
contrast, in the study described here and conducted in 
Sonora and Sinaloa, the observed maximum length of 
Pacific angel sharks was 100 cm TL, and the estimated L;, 
was 74.41 cm TL for females and 77.82 cm TL for males, 
values that correspond to 74% and 77% of the maximum 
length, respectively. 

The frequency distribution of lengths at a given time can 
vary because of different factors, including selectivity of 
fishing gear and the seasonal availability of the species in 
the area. Differences in the length-frequency distributions 
of Pacific angel sharks in the GOC have been reported. In 
this study, a peak of abundance was found for the length 
class of 85-90 cm TL for both sexes. However, for both sexes, 
Romero-Caicedo et al. (2016) reported a peak of abundance 
for individuals that were 60-90 cm TL. 

Female sharks often reach larger sizes than males 
(Castro, 2011), including sharks of the species smoothback 
angelshark (S. oculata) (Capapé et al., 2002), sawback 
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Figure 5 


The proportion of (A) female and (B) male Pacific angel 
sharks (Squatina californica) that were mature, by total 
length, among those captured from the eastern Gulf of Cal- 
ifornia in Mexico between 1998 and 2005. The dashed lines 
indicate the 95% confidence intervals of the estimated pro- 
portions from the logistic function. 


angelshark (S. aculeata) (Capapé et al., 2005), and angular 
angel shark (S. guggenheim) (Colonello et al., 2007). How- 
ever, at some point, the length structure of a stock can be 
altered by the level of historical exploitation (Stevens 
et al., 2000). Different factors, such as migration and gear 
selectivity, can influence the observed lengths in a sample. 
In particular, the effect of gear selectivity may be signifi- 
cant when interpreting life history characteristics and 
could distort growth curves and maturity ogives (Walker 
et al., 1998; Walker, 2005). 

The artisanal elasmobranch fishery in the GOC is 
intense and has been developing for several decades 
(Bizzarro et al., 2009). It is likely that the maximum sizes 
between sexes are similar in this study because Pacific 
angel sharks no longer reach the sizes they had reached 
before exploitation became extreme. The level of exploita- 
tion in the GOC does not allow individuals to reach maxi- 
mum size and affects the L5o. 

Walker (2005) suggested that reported differences in 
length at maturity may also be a result of the criteria 
adopted to define maturity. In early studies of the biology 
of sharks, length at first maturity was based mostly on the 
smallest mature individuals in the sample (Branstetter, 


1987). Such an approach fundamentally differs from the 
way L;, is estimated by using a binary logistic regres- 
sion (Walker, 2005, 2007). Technically, having insufficient 
data for the lengths of individuals that are in the transi- 
tion in maturity from immature to mature hampers max- 
imum likelihood estimation and makes estimation of the 
coefficients and 95% Cls of the logistic function challeng- 
ing (Christmann and Rousseeuw, 2001). The reproductive 
strategy of the Pacific angel shark, including vitellogenesis 
and the gestation period, can be understood only with sys- 
tematic monitoring over time (Castro, 2009). Such monitor- 
ing can be challenging because some shark species are not 
available in fishing grounds year-round. 

The incidental catch of Pacific angel sharks in artisanal 
and trawl fisheries in the GOC is concerning because of the 
low resilience of this species to fishing mortality (Cailliet 
et al., 1992). There are still gaps in our knowledge of the 
natural history of this species in the GOC; research needs 
include examination of seasonal migratory patterns and 
studies of age and growth. The estimate of L,;, from this 
study provides an important input to estimation of demo- 
graphic characteristics. The maturity ogive produced by 
using a logistic function suitably represents the develop- 
ment pattern of this species because it describes maturity 
as a gradual process. 
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Abstract—The fishery for Pacific hake 
(Merluccius productus) was established 
in the northern Gulf of California 
(NGC) in Mexico in 2000. To describe 
relevant features of the fishery and bio- 
logy of Pacific hake, we analyzed 5 years 
of data from an on board observers pro- 
gram. Using generalized additive mod- 
els, we standardized catch per unit of 
effort (CPUE), mean standard length 
(SL), and sex ratio. Also, we investigated 
by year the joint effect of maturity stage 
and sex on size distribution, biometric 
relationships, and the effect of maturity 
stage (juvenile or adult) on sex ratio. 
For all characteristics, significant inter- 
annual variability was observed. An 
average decrease was observed both in 
CPUE (52%) and mean SL (7%) during 
January—March. Adult males were 
more abundant in catch, with a male-to- 
female ratio of 1.12:1.00 (P<0.05), and 
adult females were larger than males 
(F=807.09, P<0.05). Our results support 
the previously suggested hypothesis 
that the NGC is the winter spawning 
ground for Pacific hake in the region. 
Further research is needed to deter- 
mine the status of the stock in the NGC 
and to outline specific potential man- 
agement strategies. 
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The Pacific hake (Merluccius produc- 
tus) is the most productive and eco- 
nomically important fish species on the 
west coasts of the United States and 
Canada (Hamel et al., 2015), with land- 
ings exceeding 400,000 metric tons (t) 
in recent years (Grandin et al., 2020). 
This species is distributed from the Gulf 
of Alaska to Costa Rica, including the 
Gulf of California (Lloris et al., 2003). 
In Mexico, early estimates of biomass 
of Pacific hake in the northern Gulf of 
California (NGC) ranged from 30,000 to 
65,000 t (Mathews et al., 1974; Padilla- 
Garcia, 1981), and the maximum sus- 
tainable yield was 2288-11,453 t. 
Despite the known abundance of Pacific 
hake, trawlers licensed to catch shrimp 
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and finfish species were not interested 
in commercial use of Pacific hake until 
the late 1990s (Ramirez-Rodriguez, 
2017). Total landings of Pacific hake 
have ranged from 197 t (in 2000) to a 
peak of 6628 t (in 2014), with a mean of 
approximately 2225 t for 2000-2019 in 
the NGC (Fig. 1). 

The Pacific hake, one of the most 
abundant fish species in the California 
Current Large Marine Ecosystem, plays 
a significant ecological role as a link 
between species of low (euphausiids) 
and high (sharks and seals) trophic lev- 
els (Ressler et al., 2007). Because of the 
commercial importance of this demersal 
species, its biological, ecological, migra- 
tory, and fishery traits have been widely 
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Figure 1 


Annual landings of Pacific hake (Merluccius productus) 
from the northern Gulf of California in Mexico during 2000— 
2019, in metric tons, from the Mexico National Aquaculture 
and Fishing Commission database of landing notifications. 
Data for 2019 include landings only up to 4 September. 


studied in the United States and Canada (McFarlane and 
Beamish, 1985; Smith et al., 1992; Saunders and McFar- 
lane, 1997; Benson et al., 2002; King et al., 2012). However, 
studies of Pacific hake in Mexico are few and scattered. 

The government of Mexico began issuing permits for the 
commercial catch of Pacific hake in 2018 to states around 
the Gulf of California (Sonora, Baja California, and 
Sinaloa). Currently, the only rule to control fishing effort 
is a limit of the number of permits to 80 (DOF, 2018). The 
fishing season for this species occurs, on average, during 
winter and spring (January—March) each year, and the 
largest portion of landings come from Sonora. 

Little is known about the biology and population dynam- 
ics of Pacific hake in the NGC. It has been hypothesized 
that this population belongs to a different species, the Cor- 
tez hake (M. hernandezi) (Mathews, 1985), mainly because 
of its maximum size of 107 cm total length (TL), much 
larger than that of the unfished population of what had 
been considered Panama hake (M. angustimanus), a spe- 
cies that was also called the dwarf hake, from the southern 
Gulf of California and Baja California Sur (33.5 cm TL) 
(Mathews, 1975; Balart-Pdez, 2005). Nevertheless, results 
from molecular and meristic studies indicate the presence 
of a single species, M. productus, in North America with 
different population units (Silva-Segundo et al., 2011; 
Garcia-De Leon et al., 2018). 

This study aimed to describe relevant features of the 
fishery and biology of the Pacific hake caught in the NGC. 
We report for the first time the variability of catch rates 
and population structure, by length and sex, of this species 


in this region, as well as the morphometric relationships 
(length—weight and length—length). This information can 
be used as a baseline for future stock assessments and fish- 
eries management. 


Materials and methods 
Data and sample collection 


The fishery was monitored by an on board observer pro- 
gram from 2015 through 2019 during the January—March 
fishing season with a percentage of coverage of total trips 
ranging from 5% (in 2016) to 12% (in 2019). The fleet that 
targeted Pacific hake was based in Puerto Penasco (76%) 
and Guaymas (12%), Sonora; Mazatlan, Sinaloa (8%); and 
San Felipe, Baja California (4%). During 77 trips on board 
25 commercial ships, 817 trawl tows were completed. 

The fishery used diurnal bottom trawling, and vessels 
deployed trawl net sets from 0600 to 1800. The width of 
the net mouth was 30 m on average and opened verti- 
cally (to about 4 m) as a result of a floating line in the top 
rope and a ballast chain in the lower rope of about 450 kg. 
The mesh size in the codend of nets ranged from 6.35 to 
12.70 cm (2.5—5.0 in), and the most common (89%) mesh 
sizes were 10.16 cm (58%), 7.62 cm (18%), and 8.89 cm 
(13%). The average speed and duration of trawl tows 
were 4.37 km/h (standard deviation [SD] 0.55) (2.36 kt 
[SD 0.30]) and 3.53 h (SD 1.27), respectively. The depth of 
trawl tows varied, with a range of 111-335 m and a mean 
of 253.2 m (SD 28.75). Trawl tows followed the slope of the 
Delfin Basin, always south of the Upper Gulf of California 
and Colorado River Delta Biosphere Reserve and a refuge 
for vaquitas (Phocoena sinus) in a fishing zone called the 
horseshoe by fishermen because of its shape (Fig. 2). 

On board fishing vessels, fishermen head and gut 
Pacific hake. A total of 2853 Pacific hake were weighed 
before and after processing. We estimated that an aver- 
age of 38.74% (SD 7.30) of total weight was lost during 
processing. Then, to account for the percentage of loss, the 
processed weight was multiplied by a conversion factor 
(FAO’) of 1.65 (SD 0.19). 

After processing, fish are arranged in plastic boxes that 
weigh ~70 kg, washed with seawater, and stored between 
ice layers in the cargo hold. During this process, small 
(<40 cm TL) and damaged Pacific hake are discarded. 
Nominal catch (NC) of Pacific hake per tow was estimated 
with the following equation: 


NC = CF (B x W)+D, (1) 


where CF = the conversion factor needed to estimate live 
weight from Pacific hake gutted and headed; 
B =the number of boxes stored; 
W = the weight of each box (~70 kg); and 
D = an observer's estimate of the weight of Pacific 
hake discarded. 


1 FAO. 2000. Conversion factors: landed weight to live weight. 
FAO Fish. Circ. 847, 176 p. FAO, Rome. [Available from website.] 
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Figure 2 


Map of the study area showing the fishing grounds for and catch per tow of Pacific hake (Merluc- 
cius productus) recorded in the northern Gulf of California (NGC) in Mexico during the fishing 
season from January through March in 2015-2019. The fleet that targeted Pacific hake included 
ships from Puerto Penasco, Sonora, and San Felipe, Baja California, among other localities. Trawl 
tows followed the slope of the Delfin Basin south of the Upper Gulf of California and Colorado 
River Delta Biosphere Reserve and a refuge for vaquitas (Phocoena sinus). The black solid lines 


around the Delfin Basin indicate 100-m isobaths. 


The NC was divided by the duration of the trawl tow (from 
the end of net launch until the start of net collection), to 
obtain the catch rate, or catch per unit of effort (CPUE), in 
kilograms per hour. 

Observers recorded biometrics, including standard 
length (SL) and total fresh weight, from ~60% of tows per 
trip: 98 of 131 trips in 2015, 69 of 118 trips in 2016, 93 of 
181 trips in 2017, 105 of 170 trips in 2018, and 108 of 217 
trips in 2019. They selected random samples of fish (max- 
imum of 100 fish per set) and identified them to species 
following Lloris et al. (2003). Standard length of fish was 
measured from the tip of the mouth to the caudal peduncle 
to the nearest 0.1 cm. Total fresh weight was measured 
by using a digital balance with an accuracy of 1 g. The 
Specimens were dissected to determine sex and maturity 
stages. We used a 5-stage morpho-chromatic key to deter- 
mine the maturity of Pacific hake, according to Holden 
and Raitt (1975): immature (stage 1), developing (stage 2), 
mature (stage 3), spawning (stage 4), and postspawning or 
spent (stage 5). To ease comparisons, we considered fish 
that had signs of recent, current, or imminent spawning 
(stages 3-5) to be adults. Fish assigned to stages 1 and 2 
were considered juveniles. 


Data analysis 


Catch rate standardization From all the trawl tows com- 
pleted (number of tows [n]=817), we excluded 3 tows 
because of gear malfunction, according to Maunder and 
Punt (2004), and 9 tows completed in April 2019 because 
of incomparability. A log-transformed index of catch per 
unit of effort (logCPUE) of the positive values (number 
of values=805) was used for hypothesis testing because 
these values followed a log-normal-like frequency distri- 
bution. The geometric mean was used as the index for 
mean CPUE. 

First, we developed generalized linear models (GLMs) 
to test the effect of the selected explanatory variables. 
However, we found several nonlinear effects, a result 
that violates the main assumption of GLMs. Then, 
we used generalized additive models (GAMs), which are 
semiparametric extensions of GLMs with an additive 
predictor (Guisan et al., 2002), to identify significant 
explanatory variables (speed of trawl, mesh size in the 
codend, depth, and the hour of the day, month, and year) 
and the nature of the relationships with the logCPUE. 
We considered a normal probability distribution of 
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logCPUE and used its canonical link function (identity), 
meaning that a proportional rate of change between the 
set of explanatory variables and logCPUE was expected, 
and we evaluated each model’s residuals to verify com- 
pliance with assumptions. To avoid confounding effects, 
we tested multicollinearity in explanatory variables 
through tests of correlation and analysis of variance 
inflating factors. We used only variables with variance 
inflating factors <3 (Zuur et al., 2009). Model formula- 
tions were tested according to our previous knowledge 
of the system. Stepwise addition of terms was con- 
ducted, and we kept only variables that reduced the 
Akaike information criterion (AIC) by at least 2 units 
and that had a significant increase of explained deviance 
(Marin-Enriquez et al., 2020). Also, variables with effec- 
tive degrees of freedom greater than 8 were excluded 
because they are considered highly nonlinear and dif- 
ficult to interpret (Zuur et al., 2009). The model with 
the lowest AIC was considered the best (Burnham and 
Anderson, 2004). Equal fishing power for the entire fleet 
was assumed. 


Variability of length and sex composition The length and sex 
composition of the catch of Pacific hake were assessed by 
using a 2-step approach. First, using GAMs, we repeated 
the analysis previously described to identify significant 
explanatory variables and the nature of the relationship. 
For this purpose, we used only tows for which biometrics 
were recorded (n=469) and estimated the mean SL and 
sex ratio for each one. An index of the natural logarithms 
of sex ratios (logSR) and mean SL were used as response 
variables of the GAMs. 

The second step involved a detailed analysis of the dis- 
tribution of SL in catch related to the factors of sex and 
maturity stage (juvenile or adult) by using a 2-way analy- 
sis of variance with interaction. For significant effects, we 
included the year factor to evaluate temporal consistency. 
A Tukey’s honestly significant difference test was carried 
out to identify level-specific differences. 

Also, we performed a set of chi-square tests to exam- 
ine significant departures from the null hypothesis of an 
equal sex ratio (ratio of the number of males to the num- 
ber of females: 1:1) according to Zar (1999). 


Biometric relationships The length—weight relationship 
(LWR) was estimated by fitting a power model: 


TW =a(SL)?, (2) 


where TW = the total weight (in grams); 
a = the intercept; and 
b = the slope (allometric coefficient). 


Models were fit hierarchically to each data set for Pacific 
hake caught in this study, starting with the full data 
set and continuing with data sets for each sex, maturity 
stage, and year. Parameter optimization for each model 
was conducted by minimizing the residual sum of squares 
(RSS). Statistical differences in the LWR by sex, maturity 
stage, and year were assessed by a series of “extra sum 


of squares” tests (Ritz and Streibig, 2008) defined by the 
following equation: 


_ (RSS(Bo) = RSS, (B,)) / fo - df) (3) 
RSS, / df, 


F 


where F = the Fisher’s parameter; 
df = the degrees of freedom; 
By = the nested model; and 
B, = subset of each model. 


To describe the relationships of TL and fork length (FL) to 
SL, we used linear regression models: 


y=a+b(SL), (4) 
where y = TL or FL. 


Pearson’s correlation coefficient (7) was used to evalu- 
ate the level of association of each pair of variables. The 
existence of differences between males and females for 
each relationship was determined by using analysis of 
covariance. 

All data and statistical analyses were carried out by 
using a significance level of 0.05, Microsoft Excel” 2016 
(Microsoft Corp., Redmond, WA), and statistical software R, 
vers. 3.6.1 (R Core Team, 2019). 


Results 
Catch rate standardization 


Year and the interaction of month and year were the 
most important variables (explained 10.69% of total 
deviance). Also, nonlinear effects (explained by smooth- 
ers) were found between tow speed, mesh size of trawl 
net, depth of tow, and hour of the day of tow in logCPUE 
(Table 1). The final model explained 19.9% of the total 
deviance. Levels of CPUE were higher during tows with 
speeds of 3.50—4.70 km/h (1.89-2.54 kt) and declined 
inversely with speed (Fig. 3A). Tows of trawl nets with 
mesh sizes of 7.62—8.16 cm (3.00—3.25 in) had the highest 
CPUE values, followed by tows of nets with mesh sizes 
of 9.80—10.00 cm (Fig. 3B). Values of CPUE remained 
mostly stable with depth. However, lower CPUE levels 
were observed for tows at depths >280 m, with a conse- 
quent increase in the uncertainty (Fig. 3C). The highest 
CPUE values were recorded during the morning from 
0800 to 1000 with a significant decrease beginning after 
1400 (Fig. 3D). 

Estimated mean CPUE was the highest in 2015 with 
633 kg/h (95% confidence interval [CI]: 484-830 kg/h) and 
the lowest in 2016 with 199 kg/h (95% CI: 158-250 kg/h), and 
estimates indicate less variability for 2017 (529 kg/h [95% 
CI: 368-773 kg/h]), 2018 (379 kg/h [95% CI: 310-464 kg/h]), 


? Mention of trade names or commercial companies is for identi- 
fication purposes only and does not imply endorsement by the 
National Marine Fisheries Service, NOAA. 
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Table 1 


Summary of variables in the generalized additive models used to standardize the natural 
logarithm of catch per unit of effort (logCPUE), the mean standard length (SL), and the nat- 
ural logarithm of the sex ratio (logSR) for Pacific hake (Merluccius productus) in the north- 
ern Gulf of California. The model with the lowest Akaike information criterion (AIC) was 
considered the best. Also provided are the differences in AIC values between the best model 
and each of the other models (AAICs). The explanatory variables used in the models include 
tow speed (in kilometers per hour), mesh size in the codend (in centimeters), average depth 
of trawl tow (in meters), hour of the beginning of the tow (with hours expressed as numerals 
between 1 and 24), month, year, and the interaction of month and year. Models were fit to 
data for Pacific hake caught in 2015-2019 during the fishing season from January through 
March. Asterisks (*) indicate variables that were excluded from the final models because 
their inclusion did not reduce AIC in at least 2 units along with a significant increase in the 
explained variance. The variable hour in the model for mean SL was excluded because it 
had high nonlinearity (effective degrees of freedom=8.59), considered difficult to interpret. 


Cummulative 
Response Explanatory Explained deviance 
variable variable deviance (%) explained (%) AIC 


logCPUE Null 0.00 0.00 2603.70 
Covariates 
+s(Speed) 2.43 2.43 2598.48 -5.22 
+s(Mesh) 3.69 6.12 2580.93 -17.55 
+s(Depth) 2.04 8.16 2571.97 -8.96 
+s(Hour) 1.05 9.21 2454.34 -117.63 
Factors 
+Month 1.09 10.30 2449.98 -4.36 
+Year 2.10 12.40 2439.75 -10.23 
Interactions 
+Month x Year 7.50 19.90 2380.46 -59.29 
Mean SL Null 0.00 0.00 3505.75 
Covariates 
+s(Mesh)* 0.17 3505.59 -0.15 
+s(Hour)* 4.35 3391.07 -114.68 
+s(Speed) 8.68 i 3473.78 -32.97 
+s(Depth)* 1.32 3472.30 -1.47 
Factors 
+Month 4.40 3460.61 -13.71 
+Year 4.60 3448.03 -12.53 
Interactions 
+ Month x year* 1.90 3449.56 +1.41 
Null 0.00 I 1042.23 
Covariates 
+s(Mesh)* 2.40 1042.15 
+s(Depth) 4.90 ‘ 1026.90 
+s(Speed)* 2.43 1026.00 
+s(Hour)* 0.40 985.29 
Factors 
+Month* -0.29 1030.00 
+Year 7.14 ; 1024.37 
Interactions 
+ Month x year 5.76 1013.30 


and 2019 (433 kg/h [95% CI: 321-590 kg/h]) (Fig. 4A). February (376 kg/h [95% CI: 305-454 kg/h]) and March 
In addition, a negative trend in CPUE was observed on (312 kg/h [95% CI: 253-385 kg/h]), respectively. However, 
a monthly basis (Fig. 4B). On average, January had the this pattern was not consistent across years, as indicated 
highest mean CPUE (719 kg/h [95% CI: 493-1057 kg/h]), by the relevance of the month and year interaction term 


followed by significant decreases of 48% and 52% in (Table 1). 
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Figure 3 


Partial effect plots of generalized additive models fit to the (A-D) natural logarithm of catch per unit of effort, (E) mean 
standard length, and (F) natural logarithm of sex ratio of Pacific hake (Merluccius productus) in the northern Gulf of 
California in Mexico. Explanatory variables in the models include speed of trawl tow, mesh size of trawl net, depth of 
trawl tow, and hour of the day of trawl tow (with hours expressed as numerals between 1 and 24). Data used in models 
are from Pacific hake caught in 2015-2019 during the fishing season from January through March. The gray shaded 
area in each plot represents the 95% confidence interval, and black lines just above the x-axis represent the density 
of observations for each covariate. edf=effective degrees of freedom. 


Variability of mean length 


The final model for mean SL did not include mesh size, the 
hour of the day, or depth (Table 1). This model explained 
15% of deviance and used just a smoother for tow speed 
as well as the factors of month and year without inter- 
action. The smoother for speed indicates that tows made 


slower and faster than the interval of 4.00—5.00 km/h 
(2.16—-2.66 kt) caught Pacific hake with sizes under the 
mean SL of fish caught in this study (Fig. 3E). 

The minimum annual mean SL occurred in 2015 
(45.7 cm [95% CI: 44.2—47.1 cm]), and the maximum was 
observed in 2018 (51.8 cm [95% CI: 50.5-53.2 cm]). 
Although yearly differences in mean SL were significant 
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Size distribution, sex ratio, and maturity 


The sizes of Pacific hake sampled ranged 
from 11.8 to 89.2 cm SL (15.2-108.8 cm 
TL) with a mean of 47.4 9cem SL (SD 13.7) 
(58.2 cm TL [SD 17.6]). During 2015-2017, 
the length-frequency distribution was 
bimodal. Pacific hake in the size class of 
15-25 cm SL were most common, and fish 
in the size classes of 30-45 cm SL and 
45-65-cm SL were caught in low num- 
bers. In 2018 and 2019, a more normal 
distribution of SLs was observed with 
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no clear size classes or flattened inter- 
vals. The dominant sex of Pacific hake 
at SLs >70 cm was female, and females 
had significantly larger sizes (F=807.09, 
P<0.05) than males (Table 2, Fig. 5). As 
expected, adults had significantly higher 
SL than juveniles (F=5174.49, P<0.05), 
and there was a trend across years 
(F=135.91, P<0.05). 

Adults accounted for 72% of all sam- 


T T a T lama 
2015 2016 2017 2018 2019 Jan 


Year 


Figure 4 


Observed (black dots) and predicted (black lines) geometric means for 
(A and B) catch per unit of effort (CPUE) and (C and D) mean standard 
length of Pacific hake (Merluccius productus) caught in the northern Gulf of 
California in Mexico during 2015-2019, by year and month. The gray shaded 
area in each plot represents the 95% confidence interval of the predicted mean 


from the generalized additive model. 


(P<0.05), a specific trend in mean SL was not observed 
(Fig. 4C). The monthly variability of mean SL indicates 
that the highest mean SLs occurred during January 
(49.1 cm [95% CI: 47.2-51.1 cm]) and February (50.3 cm 
[95% CI: 49.2-51.6 cm]) and that CPUE decreased 7% 
by March (47.1 cm [95% CI: 45.9—48.3 cm]), as shown in 
Figure 4D. 


Variability of sex ratio 


The final model of logSR explained 12.9% of the total devi- 
ance. The GAM used only a smoother for depth and the 
parametric effect of year as well as the interaction of month 
and year (Table 1). Results from the use of a smoother of 
depth indicate that tows with male-biased catch were con- 
ducted at depths <225 m and >275 m and that, at depths 
of 225-275 m, logSR tended to be zero (1:1 male-to-female 
ratio) or below (Fig. 3F). 

No difference was found between months in logSR, but 
a difference was observed between years. However, the 
relevance of the month and year interaction means that 
there is no consistent trend of change of logSR during the 
fishing season. 


pled Pacific hake, 70% of females, and 
73% of males. Results of the Tukey’s hon- 
est significant difference test indicate 
that all mature females were significantly 
larger than mature males, and the differ- 
ence was consistent across years. In 2015 
(adjusted P=0.99) and 2019 (adjusted 
P=0.99), no difference was detected in SL 
between sexes in juveniles (Fig. 6). 

Results from the series of chi-square 
tests indicate a male-skewed sex ratio 
for all sampled Pacific hake (male-to- 
female ratio of 1.12:1.00, P<0.05), but 
particularly for adult specimens (male-to-female ratio 
of 1.16:1.00, P<0.05). For juveniles, the sex ratio was 
1.00:1.00 (P=0.70). There was no specific annual trend 
in the sex ratio for juveniles, but for adults, there were 
significantly more males in 4 out of 5 years. In 2018, we 
observed a female-biased sex ratio for all sampled Pacific 
hake (male-to-female ratio of 0.91:1.00). However, data 
split by maturity stage indicate that the sex ratio of 
adults favored males (male-to-female ratio of 1.09:1.00) 
and that females were dominate among juveniles (male- 
to-female ratio of 0.36:1.00). Although there was no sig- 
nificant bias in the sex ratio in 2019 (P=0.09), males 
were still more abundant (Table 3). 


Biometric relationships 


We found differences in the LWR between years (F=216.55, 
P<0.05). The highest allometric coefficient (6) was esti- 
mated for 2015, followed by the lowest estimate for 2016; no 
differences were observed between 2016 and 2017 (F=1.96, 
P=0.34). For 2018 and 2019, an increasing tendency was 
observed, but differences remained present (F=86.81, 
P<0.05). Differences in the LWR by sex were found in the 
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Table 2 


Standard lengths and total weights of Pacific hake (Merluccius productus) caught in the 
northern Gulf of California during the fishing season from 2015 through 2019, by year 
and sex. Lengths and weights also are provided for all years and both sexes combined. 


n=number of samples; SD=standard d 


Standard length (cm) 


Range 


11.8-89.2 
11.8-89.1 
11.9-89.1 
11.8—79.9 
13.0—88.0 
14.0—88.0 
13.0—79.0 
14.3-89.2 
17.4-89.2 
14.3-74.1 
18.0—-88.5 
18.5-88.5 
18.0-78.1 
15.5—89.0 
16.6—-89.0 
15.5-72.0 
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eviation. 


Mean (SD) 


47.1 (13.7) 
44.8 (14.7) 
46.7 (16.1) 
43.2 (13.3) 
51.7 (12.3) 
56.4 (11.3) 
48.3 (11.8) 
48.7 (13.6) 
51.0 (14.0) 
46.6 (12.9) 
49.6 (12.0) 
52.4 (13.0) 
46.4 (10.0) 
45.8 (12.0) 
48.4 (13.9) 
43.3 (9.60) 


Total weight (g) 


Range 


22.0-9823.0 
22.0—9823.0 
25.0-9823.0 
22.0-7120.0 
26.4—7420.0 
36.0—7420.0 
26.4—4607.0 
49.0—7970.0 
79.0—7970.0 
49.0—4685.4 
47.0-8435.0 
57.0-8435.0 
47.0—4805.0 
57.2-8678.0 
66.0—8678.0 
57.2—4288.0 


10 20 30 40 50 60 70 80 90 


Standard length (cm) 


Figure 5 


Mean (SD) 


1489.9 (1016.2) 
1347.8 (1021.0) 
1558.4 (1228.1) 
1169.8 (761.5) 
1741.8 (892.9) 
2148.9 (1003.2) 
1442.5 (657.8) 
1547.4 (919.5) 
1742.0 (1036.8) 
1366.3 (751.2) 
1633.2 (1075.7) 
1931.1 (1266.4) 
1305.7 (679.6) 
1395.9 (1083.6) 
1691.5 (1358.1) 
1116.4 (616.9) 
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10 20 30 40 50 60 70 80 90 


Distribution of standard lengths of Pacific hake (Merluccius productus) caught in the northern Gulf of California in 
Mexico, by sex, during (A) 2015, (B) 2016, (C) 2017, (D) 2018, and (E) 2019. The solid and dashed lines in each panel 


represent lengths of females and males, respectively. 
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Figure 6 


Box plots of standard lengths of Pacific hake (Merluccius productus) caught in the 
northern Gulf of California in Mexico, by sex and maturity stage (juvenile [J] and 
adult [A]), during (A) 2015, (B) 2016, (C) 2017, (D) 2018, and (E) 2019. Dark and 
light gray boxes represent lengths of females and males, respectively. The thick line 
within each box indicates the median, the upper and lower parts of each box rep- 
resent the first and third quartiles (the 25th and 75th percentiles), and the whis- 
kers extending above and below each box correspond to 1.5 times the interquartile 
range. Outliers are presented as gray dots with one-tenth transparency to improve 
visualization. 


general case—with all years pooled—(F=129.16, P<0.05) 


Table 3 and consistently in every single year evaluated (Table 4; 

Fig. 7, A-E). Similarly, the model fits for 2015 and 2018 
Ratio of the number of males to the number of females, were close to each other (Fig. 7F); even so, statistical differ- 
by maturity stage Guvenile [J] or adult [A]) and by year, ences were detected (F=80.67, P<0.05). Females always had 


for Pacific hake (Merluccius productus) caught in the 
northern Gulf of California in 2015-2019. Sex ratios 
also are provided for all years and both maturity stages 


higher values of b than males. However, no significant dif- 
ferences were found in the LWR between juveniles and 


combined. The critical value for each chi-square (y”) test adults (F=0.72, P=0.97). 

was 3.81. Asterisks (*) indicate significant departures Both analyses of covariance revealed no significant effect 

from the null hypothesis of an equal sex ratio (a male-to- of sex in the TL-SL (P=0.22) and FL-SL (P=0.49) relation- 

female ratio of 1.00). M=male; F=female; and N=none. ships and a high level of correlation for both (r>0.99). The 

va ea SF Sas le ae al resulting estimates of the regression parameters for the 
Maturity Sex Observed , TL-SL relationship were 0.89 for a and 1.21 for b, and 


stage estimates of the parameters for the FL-SL relationship 


were 0.51 for a and 1.10 for bd. 


B 


ley 


Discussion 


Reports from previous studies have highlighted that the 
NGC is the winter spawning ground for the Pacific hake 
because its larvae are dominant in larval fish assemblages 
(Aceves-Medina et al., 2004; Sanchez-Velasco et al., 2009; 
Peiro-Alcantar et al., 2013). Also, the high proportion of 
adults in the catch in this zone (70%) and the sudden 
increase of catch rates in a relatively limited area found 
in this study indicate that the fishery for Pacific hake 
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Table 4 


Estimated means of the parameters a (intercept) and 6 (slope or allome- 
tric coefficient) from the regression power model used for analysis of the 
length—weight relationship of Pacific hake (Merluccius productus) caught 
in the northern Gulf of California during the fishing season in 2015-2019, 
by year, sex, and maturity stage (juvenile [J] or adult [A]). Estimates also 
are provided for all years, both sexes, and both maturity stages combined. 
Degrees of freedom (df), the coefficient of determination (r”), and the stan- 
dard error of the mean (SE) are provided for each estimate. 


Regression parameter 


Maturity 
stage df 


All 18,101 
All 8588 
9511 

5075 

1,2941 

3323 

1100 

1487 

1613 

1057 

3923 

1406 

1595 

1461 

1118 


SSSSS757 3354 


a (SE) 


0.02 (0.00) 
0.02 (0.00) 
0.05 (0.00) 
0.02 (0.00) 
0.02 (0.00) 
0.01 (0.00) 
0.03 (0.00) 
0.02 (0.00) 
0.02 (0.00) 
0.02 (0.00) 
0.01 (0.00) 
0.03 (0.00) 
0.03 (0.00) 
0.02 (0.00) 
0.02 (0.00) 


b (SE) 


2.89 (0.01) 
2.93 (0.01) 
2.63 (0.01) 
2.90 (0.01) 
2.82 (0.02) 
3.06 (0.01) 
2.78 (0.02) 
2.86 (0.02) 
2.90 (0.02) 
2.90 (0.02) 
3.02 (0.01) 
2.74 (0.02) 
2.77 (0.02) 
2.87 (0.01) 
2.90 (0.01) 


in the NGC occurs during a transient spawning aggre- 
gation (Domeier, 2012). Additionally, spawning aggre- 
gations of Pacific hake have a male-biased sex ratio, as 
has been reported for the offshore population of Pacific 
hake (Saunders and McFarlane, 1997). Several explana- 
tions have been proposed for this phenomenon in other 
hake species, including early arrival and late departure of 
males to the spawning ground, females spawning in mid- 
water and a consequently higher proportion of males at 
the bottom, and a biological strategy to ensure fecunda- 
tion (Botha, 1986; Di Giacomo et al., 1993; Pajaro et al., 
2005; Bustos et al., 2007; Korta et al., 2015). 

It is essential to note that the sampling period (and the 
fishing season for Pacific hake) begins once fishermen 
obtain high catch rates and once Pacific hake reach com- 
mercial size (>40 cm SL). The season ends once the catch 
rate and mean SL diminish to unprofitable levels. Given the 
low ex-vessel prices of Pacific hake, with an average price 
of $0.61/kg (in 2019 US. dollars) during the study period, 
the profitability of the fishery for Pacific hake depends on 
achieving high catch rates with low fishing effort (i.e., few 
fishing days). In addition, the presence of other aggregated 
and more highly valued finfish species, such as the Gulf 
croaker (Micropogonias megalops), for which the fishing 
season occurs during March—August (Ramirez-Rodriguez, 
2017; Arzola-Sotelo et al., 2018), could offer a concomitant 
explanation for the shift in the species targeted by the 
trawler fleet during the spring. 


A pattern of catch rates of Pacific hake decreasing 
during January—March has been observed in other 
studies in the Gulf of California (Mathews et al., 1974; 
Grande-Vidal, 1983; Godinez-Pérez, 2016). The reduction 
in the catch rates in this study and in these other studies 
could be related to the end of the spawning aggregation 
of Pacific hake in the NGC. The progressive absence of 
large (in length and weight) adults reduces the number of 
individuals of Pacific hake in the fishing grounds and the 
mean size of the fish caught. 

In the NGC, the high concentrations of larvae of Pacific 
hake are related to the ephemeral presence of an anticy- 
clonic eddy during winter (Sanchez-Velasco et al., 2009). 
The presence of this eddy seems to influence the spawn- 
ing areas of other fish species that inhabit the Gulf of 
California (Contreras-Catala et al., 2012). Therefore, the 
spawning areas of Pacific hake could be associated with 
the winter eddy in this zone. The influence of oceano- 
graphic structures on spawning areas has been reported 
for other hake species (Di Giacomo et al., 1993; Marrari 
et al., 2019). 

Results from previous studies indicate that at depths 
>100 m, winter temperatures (<18°C) in the NGC remain 
vertically homogeneous (Lavin and Marinone, 2003). 
Water temperature homogeneity might partially explain 
the absence of a clear depth-dependent profile of catch 
rate or mean SL because the fleet deploys trawl nets at 
depths >200 m during the daytime, when fish have been 
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Figure 7 
Standard lengths and total weights of Pacific hake (Merluccius productus) caught in the northern Gulf of California 
in Mexico during (A) 2015, (B) 2016, (C) 2017, (D) 2018, and (E) 2019. Open circles and solid lines indicate observed 
and predicted lengths and weights for females, respectively. Open triangles and dashed lines indicate observed and 
predicted lengths and weights for males, respectively. (F) Lines indicate the length—weight relationship for Pacific 
hake caught in each year evaluated. 


observed to migrate nearer to the bottom (Godinez-Pérez, 
2016). Diel vertical migration has been reported for Pacific 
hake in other regions (Hamel et al., 2015). 

Winter spawning of Pacific hake has been reported for 
other populations in the eastern Pacific Ocean (Saun- 
ders and McFarlane, 1997). However, slight differences 
exist between the spawning time and sites of Pacific 
hake from the NGC and those for Pacific hake from the 
coastal stock and from Baja California Sur. The Pacific 
hake in the coastal stock spawn during January—March 


several kilometers offshore (Agostini et al., 2006), and 
Pacific hake from Baja California Sur spawn during 
February—May, peaking in May (Balart-Paez, 2005). 
Pacific hake in NGC spawn only during December— 
May, peaking in December—February (Denton-Castillo, 
2018). The difference in spawning times and sites 
between these stocks could serve as a gene flow barrier 
that fosters the genetic isolation indicated by results 
of other studies (Iwamoto et al., 2015; Garcia-De Leon 
et al., 2018). 
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Table 5 


Transformed standard lengths (SLs) reported in the literature for Pacific hake (Merluccius productus), by source, sex, sea- 
son, region, and original type of length measurement (total length [TL], fork length [FL], or SL). Regions used in sources 
include Canada (CAN), Strait of Georgia (SG), Vancouver Island (VI), Washington (WA), Oregon (OR), California (CA), 
Baja California Sur (BCS), and the northern, central, and southern Gulf of California (NGC, CGC, and SGC). 


SL (cm) 
Range Season 


<62.3 CA 


<54.0 CAN, WA, OR, CA TL 


24.2-65.7 
9.5-61.3 
9.5-72.2 
9.5-59.5 
9.2—28.7 


Spring-summer 
Spring 


9.5-72.2 All SG 
37.7-63.6 
37.7-67.7 All VI 
37.7-59.5 
11.7-53.2 Spring 
- Winter 
7.5-33.5 All BCS 
32.5—40.7 
Summer 


Fall 
Winter 


10.4-81.5 Winter 
11.9-89.2 
11.9-89.2 


11.9-79.9 


Winter-spring 


Our results indicate that Pacific hake in the NGC pres- 
ent sex dimorphism because adult females had greater 
maximum lengths, weights, and allometric coefficients 
than males. These findings are similar to those of other 
studies that observed differences in size at age after sex- 
ual maturity (Dark, 1975; Mathews, 1975; Francis, 1983; 
Balart-Paez, 2005; King et al., 2012). Additionally, the 
maximum size that has been reported for Pacific hake was 
observed in the NGC in this study (Table 5). 

In all years except 2015, the allometric coefficient was 
lower than 3 but above the threshold of a 6 value of 2.5 
that Pauly (1984) and Froese et al. (2011) considered an 
indication of negative allometric growth. Furthermore, 
we did not find differences in the LWR between matu- 
rity stages, indicating that length and weight increase 
proportionally throughout the life of Pacific hake. In 
this regard, the variability of the allometric coefficient 
could be interpreted as seasonal growth in weight, as 
suggested by Bailey et al. (1982). They found that Pacific 
hake lose 5-10% of total weight during their spawning 
season (because of inanition) and gain 11-30% during 


Region type 


CA, OR, WA 


CGC-SGC 
15.9-88.1 All NGC 


WA, OR, CA 
WA, OR, CA 


Spring NGC 


Original 
length 
Source 


TL Best (1963) 
Alverson and Larkins (1969) 


Dark (1975) 


Mathews (1975) 
Mathews (1985) 
McFarlane and Beamish (1985) 


Beamish and McFarlane (1985) 


Stepanenko (1989) 
Saunders and McFarlane (1997) 


Balart-Paez (2005) 
Mazorra-Manzano et al. (2008) 


King et al. (2012) 


Godinez-Pérez (2016) 


This study 


their feeding season. Unfortunately, because all Pacific 
hake were sampled during winter (spawning season) and 
no stomachs were collected in our study, it was not pos- 
sible to provide evidence of seasonal changes in weight. 

Our data indicate that there was no effect of mesh size 
on mean SL. Homogeneity in mean SL can be explained 
by smaller individuals escaping because of the mesh 
size of nets and then, once the net becomes clogged with 
large individuals, selectivity of nets remaining constant. 
Likewise, the low abundance of fish in the size class of 
30-45 cm SL during 2015-2017 might be evidence of a 
size-dependent migratory pattern, ontogenetic differences 
in their spatial distribution, or interannual variability of 
abundance and distribution linked to environmental con- 
ditions (Agostini et al., 2008). 

Although these results shed light on the basic biology 
and fishery of the Pacific hake in the NGC, future studies 
should address management questions to foster sustain- 
able development. Although population characteristics, 
such as age structure, growth, maturity, and mortality, 
and fisheries ecology remain poorly understood for this 
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species in the NGC, data-limited stock assessments could 
provide temporary reference points for management. 

The maximum size of Pacific hake in the NGC (108.8 cm 
TL) is at least 3-fold of that reported for Pacific hake in the 
southern Gulf of California and Baja California Sur (33.5 cm 
TL), and information about that stock is very limited 
(Mathews, 1975; Balart-Paez, 2005; Salinas-Mayoral, 2018). 
In the short term, focusing commercial fishing for Pacific 
hake on the NGC is a potential management strategy. 

As a precautionary measure, the specification that only 
vessels licensed to take Pacific hake can take individuals 
of this species in the NGC could help to eliminate fishing 
effort over the currently allowed level that is a result of 
illegal, or incidental, catch by other ships that fish in the 
NGC. In future research, particular attention must be set 
on examining the possible adverse effects of fishing over a 
reproductive aggregation. 
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Abstract—The golden king crab (Lith- 
odes aequispinus) has_ traditionally 
been managed in the Aleutian Islands 
by using a constant catch strategy but 
interest in abundance-based manage- 
ment has emerged with the recent adop- 
tion of a size-based stock assessment 
model. Management strategy evalua- 
tion (MSE) is commonly used to quantify 
conservation and economic trade-offs of 
alternative management strategies, but 
computational constraints can impede 
the representation of all sources of 
uncertainty. We conducted a simpli- 
fied MSE for the golden king crab that 
focused on what we regarded as the 
major uncertainties, including initial 
stock abundance, future recruitment, 
estimation of mature male biomass 
and abundance, and catch implementa- 
tion error, while capturing the existing 
unique federal and state cooperative 
management framework for crab stocks 
in the Bering Sea and Aleutian Islands. 
The simplified MSE identified recruit- 
ment variability as the most important 
factor in determining overall fishery 
performance, and conservation and eco- 
nomic metrics highlight benefits of a 
15% over a 30% exploitation rate. We 
feel this simplified approach results in 
a robust analysis despite the reduced 
computational demands compared with 
those of a full MSE. Similar approaches 
can be used for other stocks, but manag- 
ers must define management objectives, 
consider stock dynamics, and identify 
factors likely to have the greatest effect 
on expected performance. 
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Strategies for calculating catch limits 
have evolved from the use of constant 
catch, constant exploitation rate, and 
fixed escapement to a focus on thresh- 
old harvest control rules (HCRs). The 
latter are considered superior because 
they include limits and targets for sizes 
of the stock and modify exploitation lev- 
els to rebuild stocks to, or sustain them 
at, healthy levels (Quinn et al., 1990; 
Zheng et al., 1993; Caddy and Mahon, 
1995; Restrepo and Powers, 1999; Punt 
et al., 2008). However, it is necessary 
to specify parameter values for HCR 
formulas (for HCRs, see Figure 1 and 
Supplementary Figure 1 [online only]) 
because each set of parameter val- 
ues will lead to different outcomes. 
Management strategy evaluation (MSE) 
is becoming an increasingly common 
tool for quantification of trade-offs 
among conservation and economic goals 
for management strategies (Punt et al., 
2016). Management strategy evaluation 
involves identifying management objec- 
tives, identifying the key uncertainties 
and representing them in a mathematical 
model of the system (hereafter referred 
to as the operating model), identifying 
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the candidate management strategies 
(the combination ofa method for estimat- 
ing stock status and an HCR), projecting 
the population forward in the operating 
model for each candidate management 
strategy, and summarizing the results 
through the use of performance metrics. 
Because MSE is computationally inten- 
sive, addressing all major uncertainties, 
including modeling, estimation, observa- 
tion, process, and implementation, can 
be impractical. As such, if MSE is to be a 
practical management tool, it must accu- 
rately represent the management frame- 
work and narrow sources of uncertainty 
down to a key subset and represent them 
in the operating model. Here, we describe 
how to employ MSE in a simplified form 
in a unique federal and state cooperative 
management framework, using the stock 
of golden king crab (Lithodes aequispi- 
nus) in the eastern Aleutian Islands as 
an example. The data-poor nature and 
complexity of the assessment of this stock 
necessitated employment of a simplified 
MSE in which the stock assessment pro- 
cess is simulated stochastically rather 
than modeled explicitly (for a similar 
approach, see Punt et al., 2008). 
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Figure 1 

Graphs of the currently implemented (A) federal and (B) state har- 
vest control rules (HCRs) for golden king crab (Lithodes aequispinus) 
in the eastern Aleutian Islands. The federal stock assessment estab- 
lishes the overfishing level on the basis of an HCR that is a function 
of instantaneous fishing mortality (F), shown in decimal form per 
year, and mature male biomass (MMB), shown in metric tons (t). The 
state HCR (HR15) determines a discrete exploitation rate, shown in 
decimal form, as a function of mature male abundance (MMA). The 
directed pot fishery closes when MMA or MMB reaches a proportion 
of 0.25 of average MMA or a proportion of 0.25 of MMB;,, which is 
35% of the unfished level of MMB (i.e., MMB given F=F,,, which is 
35% of the unfished spawning biomass per recruit), indicated by the 
vertical drops in the lines for the HCRs. The vertical drop in panel A 
does not reach zero because F attributable to bycatch can still occur in 
the groundfish fishery when the directed fishery is closed. 


The rest of this section provides a background of the 


stock of golden king crab in the eastern Aleutian Islands 
and its associated management system, which is used to 
determine which aspects are included in the operating 
model. We then outline the candidate management strate- 
gies (those currently implemented and alternative choices 
for HCRs), including some with lower and higher target 
fishing mortality rates than those previously considered 
appropriate, and compare them by using conservation and 
economic performance metrics. 


Background and fishery 
Golden king crab in the Aleutian Islands inhabit relatively 


deep water (at depths of 300-1000 m) on structurally com- 
plex habitat, such as rock and coral. They are caught by 
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using rectangular crab pots ranging in length 
from 1.2 to 3.0 m. Because of deep habitats, pots 
are deployed with longlines, with each string 
having 30-40 pots, each approximately 200 m 
apart from the next pot. Golden king crab in the 
Aleutian Islands compose what is considered one 
stock but are managed in 2 areas: east and west 
of the longitude 174°W. The fishery is male-only 
with a minimum size limit. The fishery has been 
managed by using a constant catch management 
strategy since the 1996-1997 fishing season. 
Rationalization of crab stocks in the Bering Sea 
and Aleutian Islands in the 2005-2006 fishing 
season resulted in dramatic changes in fishing 
practices, including those for golden king crab in 
the Aleutian Islands: most notably, reduced fleet 
size and increased average pot soak time (Fina, 
2005). Although the directed fishery accounts 
for most of the total fishery mortality, a small 
amount of bycatch occurs in other crab fisheries 
and the groundfish trawl and pot fisheries (<1% 
of the total number of removed animals; Leon 
et al., 2017). There was a belief by industry mem- 
bers that the stock could sustain higher levels of 
fishing intensity because no adverse effects to 
the population were detected through the use of 
fishery-dependent data under the constant catch 
management strategy (Siddeek et al., 2020). 


Assessment and management system 


A stock assessment provides inputs to HCRs 
(i.e., annual population abundance estimates 
and measures of productivity) and evaluates 
whether overfishing has occurred (i.e., whether 
total fishing mortality exceeds the correspond- 
ing overfishing level [OFL]) or whether the stock 
is in an overfished state (i.e., whether the stock 
size is below the minimum stock size threshold 
[MSST], which can be no lower than one half 
of the biomass [B] corresponding to maximum 
sustainable yield [MSY] or 0.5Bygy). The stock 
assessment for golden king crab in the Aleutian Islands 
has been based on a male-only, size-structured population 
dynamics model, which was adopted for management by 
the North Pacific Fishery Management Council in 2016. 
Unlike for other major crab stocks in the Bering Sea and 
Aleutian Islands, area-swept abundance estimates are 
not available for this stock because of the lack of a bottom 
trawl survey (NPFMC, 2018). Consequently, the assess- 
ment for the stock in the Aleutian Islands is based on 
fishery-dependent data (catch per unit of effort [CPUE] 
and catch size composition) and accounts for changes in 
fishing behavior due to crab fishery rationalization by fit- 
ting the operating model to separate sets of CPUE indices 
and by estimating fishery selectivity patterns for the pre- 
and post-rationalization periods (Siddeek et al., 2020). 
The North Pacific Fishery Management Council fishery 
management plan for the king and Tanner crabs in the Bering 
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Sea and Aleutian Islands established a state and federal 
cooperative management regime that defers crab manage- 
ment to the state of Alaska with federal oversight (NPFMC?). 
As part of the federal process, status determination criteria, 
including OF Ls and acceptable biological catches (ABCs), are 
calculated annually for crab stocks on the basis of a 5-tier 
system that accommodates varying levels of uncertainty, in 
which stocks with more biological information and greater 
assessment richness fall into lower tiers (NPFMC’). Annual 
catch levels (i.e., total allowable catches [TACs]) are deter- 
mined by the state of Alaska according to fishery regulations 
established by the Alaska Board of Fisheries, but catch levels 
and management actions need to be consistent with the pro- 
visions of the fishery management plan of the North Pacific 
Fishery Management Council, the national standards of the 
Magnuson-Stevens Conservation and Management Act, and 
other applicable federal regulations (NPFMC’). As such, it 
is necessary to consider the federal and state HCRs to ade- 
quately represent the management system. 

The federal stock assessment establishes the OFL on 
the basis of an HCR (Fig. 1A) that is a function of instanta- 
neous fishing mortality and mature male biomass (MMB), 
includes a proxy for Fyysy (F'35, the fishing mortality rate 
corresponding to 35% of the unfished spawning biomass 
per recruit; Clark, 1991, 2002) and a proxy of Bygy (B35, 
the spawning biomass corresponding to F3;). The ABC is 
then computed as 75% of the OFL. The directed fishery 
closes when MMB is <0.25Bygy. 

The state HCR (Fig. 1B, Suppl. Fig. 1B [online only]) deter- 
mines a (discrete) exploitation rate as a function of mature 
male abundance (MMA), involves a target level of average 
MMA (MMA,,,.), a threshold for opening and closing the 
directed fishery (0.25MMA,,,.), and a maximum exploita- 
tion rate on MMA (Daly et al., 2019). The outcomes from 
this HCR may be constrained by a maximum number of 
legal-sized males that are allowed to be removed, or they 
may be unconstrained. 

The final TAC is the minimum of the outcomes of the 
state HCR (when expressed as a catch) and the ABC. We 
investigated only options for the state HCR in this study 
but simulated both federal and state HCRs because, for 
example, increasing the target exploitation rate in the 
state HCR may have no effect on stock dynamics if the 
resulting catch is less than the federal ABC (i.e., is a 
lower level to prevent catch exceeding OFL). Historically, 
the outcome of the state HCR has been substantially 
lower than that of the ABC (Suppl. Fig. 2) (online only). 


Materials and methods 


The recently developed stock assessment model for golden 
king crab in the Aleutian Islands provides the input required 


' NPFMC (North Pacific Fishery Management Council). 2008. 
Final environmental assessment for Amendment 24 to the fish- 
ery management plan for Bering Sea/Aleutian Islands king and 
Tanner crabs to revise overfishing definitions, 177 p. [Available 
from website.] 


to support a state HCR that scales the target exploitation 
rate by using population abundance. However, the shift to 
abundance-based management needs to be evaluated by 
using analyses tailored to the management framework 
for golden king crab in the Aleutian Islands. We evaluated 
5 HCRs (Table 1) for the stock of gold king crab in the eastern 
portion (i.e., east of 174°W) of the Aleutian Islands by project- 
ing the population (in the operating model) forward in time 
with catches determined by the state HCRs (constrained by 
the federal ABCs). For the initial year of 2018, stock abun- 
dance by size class and operating model parameter val- 
ues were estimated by using an integrated size-structured 
assessment model based on data for 1981—2018 (Siddeek 
et al., 2020; Suppl. Materials [online only]). The 30-year pro- 
jections were replicated 1000 times. The fishing mortality 
rate for the groundfish fishery was set to the average during 
1999-2018. We chose a 30-year period for the projections 
because survival for animals in the population in 2018 would 
be <1% by 30 years, assuming an instantaneous natural 
mortality rate of 0.21/year. The trajectories of selected per- 
formance metrics appear to stabilize within ~20 years (for 
example, for trends in MMB and MMA, see Figure 2). 


Selection of uncertainties 


It is computationally impossible to consider all possible 
sources of uncertainty associated with a stock and fishery. 
Rather, the sources of uncertainty considered in analyses 
were selected because they were thought to be those most 
likely to substantially affect performance of management 
strategies (following table 3 in Punt et al., 2016): 


e Uncertainty in stock productivity and recruitment 
variability: captured with alternative values for the 
steepness (h) of the stock—recruitment relationship in 
the Ricker model (Ricker, 1954) and with the extent 
of variation and autocorrelation in recruitment in 
the stock—recruitment relationship; 

e Estimation uncertainty: captured with the extent of 
variation and autocorrelation in estimates of biomass 
and abundance and with linear and nonlinear rela- 
tionships between CPUE and stock abundance; 

e Initial stock size uncertainty: captured with alterna- 
tive specifications for initial stock abundance by size 
class and with the extent of variation and correlations 
in initial stock abundance among size classes; and 

e Implementation uncertainty: accounted for with 
the extent of variation in realized catch with regard 
to TAC. 


There are many other possible uncertainties that could 
have been but were not included in the simplified MSE, 
owing to a lack of evidence for such factors based on his- 
torical data and a lack of data with which to parameterize 
them in the operating model: 


e Process error: depensation in the stock—recruitment 
relationship and occasional catastrophic mortality 
or recruitment events; 
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Table 1 


Summary of the 5 harvest control rules (HCRs) of the U.S. government and state of Alaska evaluated for golden king crab (Lithodes 
aequispinus) in the eastern portion of the Aleutian Islands. The management strategy with zero exploitation rate (HRO) is the ref- 
erence HCR; HR10 has a maximum 10% exploitation rate with a 0.25 catch proportion cap on legal-sized male abundance, HR15 
has a maximum 15% exploitation rate with a 0.25 catch proportion cap on legal-sized male abundance, HR15U has a maximum 
15% exploitation rate without a cap on the proportion of legal-sized male abundance that can be caught, and HR30 has a maximum 
30% exploitation rate with a 0.25 catch proportion cap on legal-sized male abundance. Average mature male abundance (MMA,,,,) 
was estimated for the period from the 1985-1986 fishing season through the 2018-2019 fishing season. The HCR currently imple- 
mented by the state of Alaska for golden king crab in the eastern portion of the Aleutian Islands is HR15. F=fishing mortality rate; 


MMB=mature male biomass; OF L=overfishing level; and ABC=allowable biological catch. 


Minimum 
MMB for a 
fishery to take 


Harvest 
control 
rule 


Maximum 


Government F ABC 
Federal 


State 


Fis 0.25MMB,, 0.750FL 


e Non-stationarity: changes over time in the form and 
parameters of the stock—recruitment relationship 
and time-varying natural mortality, growth, and 
selectivity; and 

e Other factors: spatial and stock structure and time- 
varying movement. 


Key features of the operating model 


The basic dynamics are governed by this equation: 


M]y 


1,J 


Nes = De Lue ri (Cri aP Dei + Tr; ele) 


+R, (2018sts2047), (1) 


+1,j 
where N;,,; = the number of male golden king crab in 
size class at the start (1 July, the start of 
fishing year) of year t+1; 


Cti =the number of individuals retained in 
catch of the directed pot fishery for size 
class i during year tf; 


Di; =the number of individuals discarded or 
dead in the pot fishery catch for size class i 
> during year f; 
Tr,; =the number of individuals discarded or 
dead in the groundfish fishery catch for 
size class i during year f; 
X;; =the probability of animals in size class 7 
growing into size class j during the year; 
y, = the time from 1 July to the midpoint of the 
fishery period during year ¢; 
M = the instantaneous rate of natural mortality 
(assumed to be 0.21/year for all size classes 
and over time; Siddeek et al., 2020); and 


Maximum 
exploitation 
rate 


Minimum 
MMA for a 
fishery to take 


Catch proportion 
cap on legal-sized 
male abundance 


Catch 
limit 


0.25MMA,,,. 
0.25MMA,,,,. 
0.25MMA,,,,. 
0.25MMA,,. 
0.25MMA,,. 


Ri; =the recruitment to size class j during 
year t+1. 


The equations used to compute future retained and dis- 
carded catches are provided in Equations Al—A4 in 
Supplementary Materials (online only). 

Future recruitment is an essential part of the operating 
model. However, fitted stock—recruitment models often fail 
to show the link between spawning individuals and recruits 
for many species (Subbey et al., 2014), including for most 
commercially important crab stocks in the Bering Sea and 
Aleutian Islands. Biomass of mature females is a general 
choice for an index of egg production for finfish species (e.g., 
Martell et al., 2008; Punt et al., 2008; Subbey et al., 2014). 
However, MMB has been adopted as a proxy for reproduc- 
tive output for stocks of golden king crab in the Aleutian 
Islands, largely as a result of uncertainties related to iden- 
tifying the component of the mature population that par- 
ticipates in mating and to defining optimal sex ratios and 
because fisheries have a male-only retention requirement 
(NPFMC’). The MMB is computed by using this equation: 


MMB, = dete ore 


=(Cy,5+Dij+ Tr, pe hy, (2) 
where y’ = the time from 1 July to 15 February in the fol- 
lowing year (NPFMC’); 
w; = the weight for size class j (Siddeek et al., 2020); 
m=the lowest size class with mature animals 
(i.e., maturity is assumed to be a knife-edged 
function of size at the carapace length [CL] of 
111 mm; Daly et al., 2019); and 


n = number of size classes. 


384 


20,000 


15,000 


10,000 


Year 


Fishery Bulletin 118(4) 


MMA (millions) 
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HCR — HRO --- HR10 --- HR15 -- HR15U --- HR30 


Figure 2 


Median mature male biomass (MMB) and median mature male abundance (MMA) of golden king crab 
(Lithodes aequispinus) in the Aleutian Islands under 5 harvest control rules (HCRs), with the initial state of 
the stock set at (A and B) healthy or (C and D) overfished, for scenario 1 of the operating model in which a 
linear relationship between catch per unit of effort and selected abundance is assumed. Metric tons (t) is the 
unit for MMB, and number of crab is the unit for MMA. Values are based on a 30-year projection period that 
begins with 2018. The stock is projected from 2 initial levels of abundance: a healthy state (i.e., MMB.;¢/ 
MMB;;=1.55, where MMBoo 3 is MMB in 2018 and MMB,, is 35% of the unfished level of MMB) and an over- 
fished state (i.e., MMB..1:/MMB,;=0.50). Horizontal dashed lines indicate the MMB,, and average MMA 
(MMA,,,.) thresholds. The HCRs include HRO, the reference HCR with a zero exploitation rate; HR10, with a 
maximum 10% exploitation rate and a 0.25 catch proportion cap on legal-sized male abundance; HR15, with 
a maximum 15% exploitation rate and a 0.25 catch proportion cap on legal-sized male abundance; HR15U, 
with a maximum 15% exploitation rate and no cap on the proportion of legal-sized male abundance that can 
be caught; and HR30, with a maximum 30% exploitation rate and a 0.25 catch proportion cap on legal-sized 


male abundance. Data used in the model are for golden king crab in 1981-2018. 


Mature male abundance is determined from Equation 2 
without w;. 


Parameterization of the operating model 


The population dynamics model has seventeen 5-mm size 
classes over the size range of 101-185 mm CL. The last size 
class (181-185 mm CL) is a plus group that includes all golden 
king crab larger than 185 mm CL. The values for the param- 
eters of the model were estimated on the basis of fitting the 
model to available data (Siddeek et al., 2020; for some of the 
estimated parameters, see Supplementary Figure 3 [online 
only]). Although Bayesian or Monte Carlo methods could have 
been used to estimate distributions for the parameters of the 


operating model, these methods were not pursued because 
they were found to be computationally infeasible. The param- 
eters of the operating model differed depending on whether 
CPUE was assumed to be proportional to selected abundance 
or the square root of selected abundance because the oper- 
ating model was fitted to the available monitoring data for 
golden king crab in the Aleutian Islands (for the results for 
linear relationship [hereafter referred to as the linear choice], 
see Tables 2 and 3 and Supplementary Tables 1 and 2 [online 
only]; for the results for the nonlinear relationship [hereaf- 
ter referred to as the nonlinear choice], see Supplementary 
Tables 3-6 [online only]). Tables and figures for the linear 
choice are provided (selected tables and figures for the nonlin- 
ear choice are given in Supplementary Materials [online only]). 
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The value of h was determined through the use of a 
2-step procedure. The first step was to calculate F;,, by using 
an analysis of spawning biomass per recruit (for a plot of 
spawning potential ratio [SPR] versus fishing mortality [F], 
where the SPR for a given level of fishing mortality is the 
ratio of the spawning biomass per recruit at that level of 
fishing mortality to the spawning biomass per recruit in an 
unfished state, see Supplementary Figure 4A [online only]). 
The value of F is related to the total fishing mortality for 
the directed pot fishery when computing SPR, with fishing 


Table 2 


The factors considered during the evaluation of the perfor- 
mance of harvest control rules for golden king crab (Lith- 
odes aequispinus) in the eastern portion of the Aleutian 
Islands. The 3 values are lower limits; the model estimates 
for 2018, indicated with asterisks (*); and upper limits con- 
sidered in the evaluation of the performance of harvest con- 
trol rules. The values are based on the operating model in 
which a linear relationship between catch per unit of effort 
and selected abundance is assumed. Data used in the model 
are for golden king crab in 1981-2018. MMB=mature male 
biomass; MMA=mature male abundance. 
Factors Values 

0.600, 0.729%, 1.200 
0.100, 0.274*, 0.800 
0.000, 0.455*, 0.900 
0.000, 0.039%, 0.100 


Steepness (h) 

Recruitment variation (Gp) 

Autocorrelation in recruitment (pp) 

Extent of annual catch 
implementation error (6) 

Extent of MMB estimation error (6g) 0.000, 0.214*, 0.300 

Autocorrelation in MMB estimation 0.000, 0.700, 0.900 
error (pp) 

Extent of MMA estimation error (6x) 0.000, 0.216*, 0.300 

Autocorrelation in MMA estimation 0.000, 0.700, 0.900 
error (Py) 


mortality attributable to bycatch in the groundfish fishery 
set to the estimated average during 1999-2018 (although 
this source of mortality is negligible). The second step was 
to calculate MMB for various h values by projecting the pop- 
ulation model forward deterministically given F=F, and to 
find h such that MMB=0.35xMMB, (where MMB, is the vir- 
gin mature male biomass; Suppl. Fig. 4B [online only]). 

The selected value of h ensures that fishing at F;,, pro- 
duces MMB,,, which is equivalent to 0.35xMMB, (a proxy 
for Bysy), and a proxy MSY (see Supplementary Figure 5 
[online only]). Figure 3 shows the fitted (deterministic) 
Ricker stock—recruitment relationship when h was set 
by using this procedure (h=0.729, Rj=2.528 million crab, 
MMB,=18,862 metric tons) along with the MMB,_, (MMB 
estimated 8 years prior to year ¢t) and recruit (R,) pairs 
of data points from the stock assessment model. Figure 4 
depicts the estimated number of recruits from the model, 
along with a 95% confidence interval, indicating that 
uncertainty in these estimates increased over time. 


Harvest control rules 


Federal harvest control rule The federal HCR is based on 
MMB and is needed to determine the OFL and the ABC for 
each year of the projection period. The OFL fishing mortal- 
ity (For) is determined by using this equation (NPFMC’): 


Fs, if MMB,>MMB,,, 

MMB, nee 

MMB;, 
l-a 


Fa =y Fa if 0.25MMB,.<MMB,<MMB,,, and (3) 


0 if MMB,«0.25MMBy,, 


where o = 0.1 (a preset value). 


The OFL of total catch is calculated by applying Foy, to 
male stock abundance through the use of Equations A1—A3 
in Supplementary Materials (online only) and by summing 


Table 3 


The reference points used in the evaluation of candidate harvest control rules for golden king crab 
(Lithodes aequispinus) from the eastern portion of the Aleutian Islands: a proxy for biomass corre- 
sponding to maximum sustainable yield (35% of the unfished level of mature male biomass [MMB,,]), 
the fishing mortality rate corresponding to 35% of the unfished spawning biomass per recruit (F';), 


mean catch, average mature male abundance (MMA 


ave 


), and mean catch per unit of effort (CPUE). 


The estimates are based on the operating model in which a linear relationship between CPUE and 
selected abundance is assumed. Data used in the model are for golden king crab in 1981-2018. 


Reference point Estimate 


MMB;; 
P55, 0.656/year 


6601.61 t 


Mean catch 1492.82 t 


Basis of estimation 


Assessment model 
Assessment model 
Period from the 2005-2006 fishing season through the 


2018-2019 fishing season (post-rationalization period) 


MMA 5.45 million crab 


ave 


Period from the 1985-1986 fishing season through the 


2018-2019 fishing season (estimation period) 


Mean CPUE 31.70 crab/pot lift 


Period from the 2005-2006 fishing season through the 


2018-2019 fishing season (post-rationalization period) 
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The currently adopted state HCR is HR15, 
with a maximum exploitation rate of 15% and a 
0.25 catch proportion cap on abundance (num- 
ber) of legal-sized males. An alternative HCR is 
HR10, with a maximum exploitation rate of 10% 
and a 0.25 catch proportion cap on legal-sized 
male abundance; HR15U is an alternative HCR 
with a maximum exploitation rate of 15% with- 
out any cap on the proportion of legal-sized male 
abundance that can be caught, HR30 is an alter- 
native HCR with a maximum exploitation rate of 
30% and a 0.25 catch proportion cap on legal-sized 
male abundance, and HRO is a reference HCR 
with no directed fishery (i.e., zero exploitation 


5 + Steepness=0.729 


MMB35=6602 | 


No. of recruits (millions) 


MMA aye=5-45 


5000 rate) (Table 1, Fig. 1, Suppl. Fig. 1 [online only]). 
MMB (t) The state HCR is based on MMA and computes 
the catch proportion for year t, HR,, as follows: 
Figure 3 
The stock—recruitment relationship from the Ricker model fitted to the : 
stock and recruitment data for golden king crab (Lithodes aequispinus) i eta ae 
in the eastern Aleutian Islands during 1981-2018. The vertical dashed HR, = | sills | if 0.25MMA,,,<MMA,<MMaA,,,, and (6) 
lines indicate the respective reference points: 35% of the unfished level MMA.y< 
of mature male biomass (MMB,,;), which is a proxy for biomass cor- 0 if MMA,«<0.25MMA,,., 


responding to maximum sustainable yield, and average mature male 
abundance (MMA,,,,). The steepness value is the estimate for 2018. In 
the operating model used in this analysis, catch per unit of effort is 
assumed to be proportional to selected abundance. t=metric tons. 


where x = the maximum exploitation rate; and 
MMaA,,,. = the average MMA during 1985-2018. 


ave 


The exploitation rate from Equation 6 is con- 
verted into an instantaneous fishing mortality 
over the predicted catches in weight by source of mortal- rate (F, in year t) by solving this equation: 


ity, in other words, by using this equation: 


m f{~, a ae F. x ge x 
OFL, =¥i(c jj + D', jwjt+Tr ol03)> (4) HR, == =e * (7) 
t 
where the prime symbols indicate that the catch 
and bycatch are calculated by using F=F'op, 4. 
The ABC (which applies to all sources of fishery- 
related mortality) is calculated (NPFMC, 2018) 


with this equation: 


ABC, = 0.75 x OF L,. (5) 


SY 
a 


State harvest control rule options We compared 
5 candidate state HCRs, with the aim to main- 
tain consistency with state of Alaska commercial 
fishery regulations, the Alaska Board of Fish- 
eries policy on king and Tanner crab resource 
management (ABF”), the North Pacific Fishery 
Management Council fishery management plan 
(NPFMC!), and National Standards 1 and 2 of 
the Magnuson-Stevens Fishery Conservation 
and Management Act (NPFMC!). The candidate 
HCRs were informed by catch policies for other 
stocks of king crab species in the Bering Sea and 


No. of recruits (millions) 


1990 2000 2010 
Year 


Aleutian Islands, historical exploitation rate esti- 
mates for the fishery for golden king crab in the 
Aleutian Islands, and stakeholder input. 


? ABF (Alaska Board of Fisheries). 1990. Policy on king 
and Tanner crab resource management. Policy no. 90- 
04-FB. [Available from website.] 


Figure 4 


The number of recruits of golden king crab (Lithodes aequispinus) 
in the eastern Aleutian Islands for 1981-2019 estimated with the 
stock assessment model. The gray shaded area indicates the 95% con- 
fidence interval. In the operating model used in this analysis, catch 
per unit of effort is assumed to be proportional to selected abundance. 
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where S" = the total selectivity (a curve for selection of 
all susceptible sizes of retained and discarded 
crab to pot gear, fixed at the assessment model 
estimates); and 
Z' =the instantaneous total mortality during 
year t (including components for the directed 
pot fishery and bycatch in the groundfish 
fishery). 


A cap on the proportion of legal-sized male abundance 
that can be caught is included in the state HCRs for sev- 
eral Alaska crab stocks, such as the stock of red king crab 
(Paralithodes camtschaticus) in Bristol Bay (Zheng et al., 
1995). The HR, may therefore be further modified for con- 
servation purposes by constraining the predicted catch (in 
numbers) by the directed pot fishery to not exceed 25% of 
the legal-sized male abundance (crab size >136 mm CL; 
Siddeek et al., 2020), in other words, 


yeg@s js0.255" Ny, 


where L = the legal minimum size of golden king crab in 
the Aleutian Islands. 


The predicted catch (in weight) by the directed pot fishery 
cannot exceed the retained component of the ABC (cur- 
rently adopted precautionary state management policy is 
to avoid exceeding the federal fishery management limit, 
the OF L), in other words, 


ers A 
Yale < 0.75 x Deal t, jYj- 
Representing uncertainties 


Uncertainty about future recruitment Many types of stock— 
recruitment relationships can be fitted to the results 
from an assessment, but the Ricker model (Ricker, 1954) 
was chosen for this study given its previous use for king 
crab species (Lithodidae spp.) (e.g., Zheng et al., 1995; 
Bechtol and Kruse, 2009). Future recruitment is gener- 
ated with variation and temporal autocorrelation, with 
this equation: 

MMB, x -125n(5n) Ska o,-oR 
15 Fo wap, ° 

0 


€, =Pretatyl- Pre» and 


e~ N (0.0%), 


R, 


where R, = the number of recruits at unfished equilibrium; 
MMB, ,, = the MMB (in metric tons) in year t-k given 
a k-year lag between spawning and recruits 
entering the model (k=8; Daly et al., 2019); 
MMB, = the unfished MMB; 
h = the steepness parameter; 
Pp =the extent of autocorrelation in the recruit- 
ment deviations; 
Op = the standard deviation of recruitment; and 


Q; = a normalized gamma function that determines 
the distribution of recruits to each size class: 


1,+2.5 
ii gamma(x| o.,, B,.)dx 
1,-2.5 
Oe and (9) 
J n 1}+2.5 A 
Dealers gamma(s| a, B,)dx 
ax 
a,—1 By 
gamma(x |o,, B,) = x, 
B, NG.) 


where @, = a parameter of the gamma distribution; 
B,. = a parameter of the gamma distribution; 
1; = the midpoint of size class j; and 
n = the number of recruiting size classes, fixed to 5. 


Uncertainty about the state of the stock in 2018 For simplic- 
ity, most parameter values were fixed at their best esti- 
mates. However, uncertainty was introduced to the size 
composition for the first projection year (2018), according 
to the following equation: 


oe 8; ze 
Noo18,; =N2018,je ~ ? (10) 
hg Ed MN (0, YW). 


where N 2018, ; = the estimate of the number of males in 
size class j at the start of 2018; 
69018; = the standard error of the logarithm of the 


estimate of N 2018, ;; and 


V =the variance—covariance matrix for the 
numbers by size class at the start of 2018. 


Estimates of N201s,; and V were obtained from the assess- 
ment model. Stock status in each projection was deter- 
mined by using the best estimates from the assessment 
model rather than by sampling parameter vectors from a 
posterior or a bootstrap distribution, as would commonly 
be done in a full MSE. 


Uncertainty when applying the federal and state harvest 
control rules Uncertainty in the estimates of MMB and 
MMA are accounted for by replacing MMB and MMA in 
Equations 3 and 6 as follows: 


95 
uD 


: 5 
MMB¢™"e4 — MMB,e 2, (11) 
5: = Pp. t+ yl - Pee, and 
D, ~ N(0, on): and 


oN 


MMAS*™=ted — MMAe * 2 , (12) 
MO, =PyO + yl - Pn St, and 


DO, ~ N(0, ox), 


where MMBetimated = the estimate of MMB for year f; 
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MMA‘“™s'ed = the estimate of MMA for year f; 
MMB, and MMA, = the true MMB and MMA during 
year ¢ in the operating model; 
Pp and py = the extents of autocorrelation in 
stock status estimation error; and 
Op and On = the extents of estimation error. 


The values for 6g and oy are set to the standard devia- 
tions of the logarithms of the estimates of MMB and MMA 
for 2018 rather than those of MMB/MMB,, and MMA/ 
MMA,,,, because the estimates of average MMB and MMA 
are precise. The autocorrelation in estimation error cannot 
be obtained from the assessment model; therefore, a range 
of plausible values are considered in the analyses. 


Implementation error 


The fishery does not catch the TAC exactly; therefore, 
implementation error is introduced as follows: 


Crewe Caireaeyand (13) 


2 
Tt ot N(0, 06,5) 


where (Ofer = the true catch for animals in size class j 
during year f; 
C;,;= the expected catch of animals in size class 
j during year t based on the simulated fish- 
ing mortality from the HCR; and 
Oc ¢,; =the standard deviation of the differences 
between TACs and actual catches for size 
class 7 based on the standard deviation of 
the differences between the TAC and total 
landed catches, 6c, in other words, 


SOC ia aa ae (14) 
yee t,j 
where ye Ct,j =the sum of the expected retained catches 
a of all size classes during year t. 


9G; 


Simulation design 


The design of the simplified MSE involved combining 
levels for each of the uncertainties. In total, 53 scenarios 
based on the selected uncertainties were considered for 
each relationship between CPUE and selected abundance 
(Table 2; Suppl. Tables 1—5 [online only]). 

Scenario 1 was based on the best estimates of the param- 
eters, and the specifications of this scenario are indicated by 
asterisks in Table 2 and Supplementary Table 5 (online only). 
Scenario 1 was also based on the middle level of autocorrela- 
tion in error when MMB and MMA were estimated (for the 
full list of scenarios and specifications, see Supplementary 
Tables 1 and 2 [linear choice] [online only] and Supplementary 
Tables 3 and 4 [nonlinear choice] [online only]). Scenarios 
2-17 involved changing the value of one of the parameters 
of scenario 1, and scenarios 18-53 changed the value of more 
than one parameter. The scenarios did not explore all possible 


combinations of parameters owing to computational and pre- 
sentational limitations. 

Two options (applied separately for each of the linear and 
nonlinear choices) were considered for the size structure at 
the start of the projection period (1 July 2018): 1) estimate 
in the assessment model (i.e., MMB/MMB,5=1.55; Siddeek 
et al., 2020) and 2) the MSST (0.5MMB,,). The second option 
was implemented by increasing the fishing mortality rate 
on the size structure for 2018 such that MMB approached 
0.5MMB;3;. Performances of candidate HCRs were eval- 
uated by projecting the stock from the initial abundance 
levels from these 2 options, in other words, a healthy state 
(MMB>MMB;;) and an overfished state (MMB=0.5MMB,,;). 


Performance metrics 


We considered conservation and economic criteria when 
evaluating the candidate HCRs. The conservation criteria 
were 1) the probability (across simulations and the entire 30- 
year period) of the stock being below MSST (..e., a threshold 
for being overfished), 2) the probability of total catch being 
greater than OFL (i.e., a threshold for overfishing occur- 
ring), 3) the probability of total catch being greater than 
ABC, and 4) the probability that MMB is less than MMB,,. 
The economic criteria were 1) the probability of fishery clo- 
sure, 2) the average annual catch (across simulations) for 
the directed fishery, 3) the annual variability of catch in the 
directed fishery, 4) the probability of retained catch being 
less than mean retained catch during 2005-2018 (post- 
rationalization period), 5) the average CPUE, 6) the prob- 
ability of CPUE being less than mean CPUE for the period 
2005-2018, 7) fishing effort (number of pot lifts, as approx- 
imated by using Catch/[CPUEx0.00195], where mean crab 
weight is assumed to be 0.00195 metric tons on the basis 
of unpublished data [Ben®] on the retained catch in 2018), 
and 8) the probability of MMA being less than MMA,,, (an 
average value for the period 1985-2018), an indication of 
whether the exploitation rate for a given management 
strategy reaches the maximum allowable exploitation rate 
in expectation. The reference points for comparison with 
simulation results are listed in Table 3 (linear choice) and 
Supplementary Table 6 (nonlinear choice) (online only). 

In addition, the time series of MMB, MMA, catch, effort, 
and catch variation were summarized as follows: 


e Time trends in median (over simulations) MMB, 
MMA, catch, and effort over the projection period; 

e Rebuilding time from the _ overfished level 
(0.5MMB,,;) to MMB;, and the rebuilding time from 
the corresponding MMA to MMA,,,; and 

e Median (over simulations) annual catch variability 
(AAR) computed with this equation: 


WAR] lesen) (15) 
XC 


3 Ben, D. 2018. Unpubl. data. Div. Commer. Fish., Alaska Dep. 
Fish Game, 351 Research Ct., Kodiak, AK 99615. 
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The 5 HCRs were ranked within each conservation and 
economic performance criteria leading to a risk matrix 
that involved grouping these performance metrics into 3 
categories: conservation, catch, and catch stability. Overall 
ranks within each category were based on average ranks 
across all criteria. 


Results 
Performance of harvest control rules (scenario 1) 


The time trajectories of MMB and MMA for scenario 1 
(i.e., the best estimates of the parameters when CPUE 
was assumed to be linearly proportional to selected abun- 
dance; Suppl. Table 1 [online only]) stabilized toward the 
end of the 30-year projection period in median terms for 
the 5 HCRs. As expected, HR10, HR15, HR15U, and HR30 
led to much lower MMB and MMA than HRO. The alter- 
native policy of HR30 resulted in slightly lower MMB (Fig. 
2, A and C) and MMA (Fig. 2, B and D) than HR10, HR15, 


20 30 


HCR — HR10 


Effort (no. of pot lifts) 


Q 
= 
= 
fe) 
Q 
— 
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fe} 
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= 
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and HR15U, with the latter 2 HCRs having very similar 
trajectories of MMB and MMA. The MMB and MMA tra- 
jectories remained above MMB,, and MMA,,,, respec- 
tively, throughout the projection period when the 
simulation started above MMB,, and reached the respec- 
tive reference points toward the end of the projection 
period when the simulation started in an overfished state 
(Fig. 2). The estimates of MMB and MMA are qualita- 
tively similar, a result that was expected given MMB is 
MMA multiplied by average weight, which did not change 
much over time. The results are similar for scenario 1 
when CPUE was assumed to be proportional to the square 
root of abundance (Suppl. Tables 7—9) (online only). 

Catch and effort stabilized (in median terms) toward 
the end of the 30-year projection period for HR10, HR15, 
HR15U, and HR30. However, the time trajectories of 
catch and effort differed depending on the initial state. 
The median catch declined over time when MMB, 1, was 
greater than MMB,, (Fig. 5A) but increased over time for 
when MMB, 9,3 was equal to 0.5MMB,, (Fig. 5C). Effort 
was almost constant over the projection period when 


20,000 


~- HR15 --- HR15U -- HR30 


Figure 5 


Median catch and median effort of the directed pot fishery for golden king crab (Lithodes aequispinus) in the 
eastern Aleutian Islands under 4 harvest control rules (HR10, HR15, HR15U, and HR30), with the initial 
state of the stock set at (A and B) healthy or (C and D) overfished, for scenario 1 of the operating model in 
which a linear relationship between catch per unit of effort and selected abundance is assumed. Values are 
based on a 30-year projection period that begins with 2018. The stock is projected from 2 initial levels of 
abundance, measured in mature male biomass (MMB): a healthy state (i.e., MMBo,./MMB,,=1.55, where 
MMB; is MMB in 2018 and MMB,, is 35% of the unfished level of MMB) and an overfished state (i.e., 
MMB 4;:/MMB;;=0.50). Data used in the model are for golden king crab in 1981-2018. For details about the 


harvest control rules, see Table 1. t=metric tons. 


390 


Fishery Bulletin 118(4) 


MMB. 1, Was greater than MMB,, (Fig. 5B) but increased 
over time when MMB, .,;3, was equal to 0.5MMB,,; 
(Fig. 5D). The alternative of HR30 led to slightly higher 
catches but substantially higher effort when MMB, ;. 
was greater than MMB.,,. The alternative of HR10 led to 
lower catches than HR15, HR15U, and HR80. The alter- 
native of HR15U led to negligibly higher catches than 
HR15 (Fig. 5A). Results were similar for scenario 1 when 
CPUE was assumed to be proportional to the square root 
of abundance (Suppl. Table 8) (online only). 

The probabilities of the stock of golden king crab in 
the Aleutian Islands being overfished (MMB<MSST) and 
severely overfished (VUMB<0.5MSST) and of overfishing 
occurring (i.e., catch exceeding OFL) were zero, and the 
probability of this stock being below MMB,; was <0.02 
for all policies when MMB, ;, was greater than MMB;,; 
when CPUE was assumed to be linearly proportional 
to abundance (Table 4). The probability of this stock of 
golden crab being below MMB,, was higher for HR15, 
HR15U, and HR30 when CPUE was proportional to the 
square root of selected abundance, but the trends were 
similar to those with the linear relationship of CPUE to 
abundance (Suppl. Table 7) (online only). The fact that the 
probability of overfishing occurring was zero across poli- 
cies was expected, given the constraint that the predicted 


catch (in weight) in the directed fishery could not exceed 
the retained catch component of the ABC. Interestingly, 
the probabilities of the stock of golden king crab in the 
Aleutian Islands being severely overfished and of over- 
fishing occurring were also zero when the stock was ini- 
tially overfished. Probabilities of this stock being below 
MMB.,,, were greater for all policies when it was initially 
in an overfished state, yet probabilities of this stock stay- 
ing overfished during the projection period were <0.031 
for all policies, indicating a resiliency of this population 
in the simulations. 

Results from consideration of economic criteria for 
HRO were largely moot; therefore, only HR10, HR15, 
HR15U, and HR380 were ranked for these metrics. The 
alternative of HR10 resulted in lower catch in compari- 
son with that of the rest of the HCRs but led to improved 
performance for other economic metrics, such as reduced 
effort (i.e., fewer pot lifts needed to achieve a low TAC; 
see Tables 4 and 5 for the linear choice and Supplemen- 
tary Tables 7 and 8 [online only] for the nonlinear choice). 
Most economic criteria were similar between HR15 and 
HR15U for both initial conditions, indicating that the 
catch limit on legal-sized male abundance had little 
overall effect (see Table 5 for the linear choice and Sup- 
plementary Table 8 [online only] for the nonlinear choice). 


Table 4 


Conservation performance metrics for golden king crab (Lithodes aequispinus) in the Aleutian Islands, 
with the initial state of the stock set at healthy or overfished, for scenario 1 (based on best parameter esti- 
mates) of the operating model in which a linear relationship between catch per unit of effort and selected 
abundance is assumed. Values for the 5 harvest control rules (HCRs) evaluated by conservation criteria, 
HRO, HR10, HR15, HR15U, and HR30, are probabilities that the estimated quantity, such as mature male 
biomass (MMB) or total catch, is above or below the associated reference point, such as minimum stock 
size threshold (MSST), overfishing level (OFL), allowable biological catch (ABC), or 35% of the unfished 
level of MMB (MMB,,;), calculated for the last 10 years of the 30-year projection period, which begins with 
2018. For example, values for MMB<MSST are the probabilities that MMB is below MSST. Harvest control 
rules were ranked among each other for each performance metric (ranks are given in parentheses; ranks 
are the same for HCRs if probabilities are the same for those HCRs). Total catch is the catch retained in 
the directed pot fishery plus the discard mortality in the directed fishery and the bycatch mortality in the 
groundfish fishery. The stock is projected from 2 initial levels of abundance, measured in MMB: a healthy 
state (i.e., MMBoo13/MMB,,=1.55, where MMB.o1, is MMB in 2018) and an overfished state (i.e., MMBo91¢/ 
MMB,;=0.50). For details about the HCRs, see Table 1. Byjgy=the biomass corresponding to maximum sus- 


tainable yield. 
Metric 


Healthy 


Description HRO 


Overfished 
Severely overfished 
Overfishing (OFL) 
Overfishing (ABC) 
Below Bysy 

Overfished 
Overfished 
Severely overfished 
Overfishing (OFL) 
Overfishing (ABC) 
Below Bysy 


MMB<MSST 
MMB<0.5MSST 
Total catch>OFL 
Total catch>ABC 
MMB<MMB,, 


MMB<MSST 
MMB<0.5MSST 
Total catch>OFL 
Total catch>ABC 
MMB<MMB;, 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 


HR10 HR15 HR15U HR30 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.001 (3) 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.001 (3) 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.016 (5) 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.063 (4) 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.141 (5) 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 


0.000 (1) 
0.000 (1) 
0.000 (1) 
0.000 (1) 
0.063 (3) 
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Table 5 


Economic performance metrics for golden king crab (Lithodes aequispinus) in the Aleutian Islands, with the initial 
state of the stock set at healthy or overfished, for scenario 1 (based on best parameter estimates) of the operating 
model in which a linear relationship between catch per unit of effort and selected abundance is assumed. Values 
for the 4 harvest control rules (HCRs) evaluated by economic criteria, HR10, HR15, HR15U, and HR30, are prob- 
abilities that the estimated quantity is above or below the associated reference point (e.g., MMA<0.25MMA,,,, 
indicates the probability that mature male abundance [MMA] is below 25% of average MMA [MMA,,,,]), calculated 
for the last 10 years of the 30-year projection period, which begins with 2018. The exceptions are for values of catch 
(mean in metric tons), CPUE, (given as the number of crab per pot lift), and effort (given as the number of pot 
lifts). Harvest control rules were ranked among each other for each performance metric (ranks are given in paren- 
theses; ranks are the same for HCRs if probabilities are the same for those HCRs). Catch is the catch retained in 
the directed pot fishery. The stock is projected from 2 initial levels of abundance, measured in mature male bio- 
mass (MMB): a healthy state (i.e., MMBoo19/MMB,;=1.55, where MMB, o;, is MMB in 2018 and MMB,, is 35% of 
the unfished level of MMB) and an overfished state (i.e., MMBo9)3/MMB,,=0.50). For details about the HCRs, see 
Table 1. Catcha,.7i,¢=historical average catch; CPUE,,,.4;,.=historical average CPUE; and MMA,,,4;,,=historical 


average MMA. 


Metric Description or unit 


HR10 HR15 HR15U HR30 


Healthy 
Fishery closure 


MMA<0.25MMA,,,. 
Catch Mean 


Catch variability 
Relative catch Catch<Catcha,.nist 
Mean number 
CPUE<CPUE gy crist 
Number of pot lifts 
Stock status MMA<MMA,, Hist 
Overfished 
Fishery closure MMA<0.25MMA,,,. 
Catch Mean 
Catch variability 
Relative catch 
CPUE, Mean number 
CPUE, CPUE<CPUE yerist 
Effort Number of pot lifts 


Stock status MMA<MMA yy crist 


Catch<Catcha, nist 


The alternative of HR30 was the most aggressive policy, 
but with only marginally higher catches than HR15 and 
HR15U, and HR15U had the most variable catches (see 
Table 5 for the linear choice and Supplementary Table 8 
[online only] for the nonlinear choice). 


Performance of harvest control rules (all scenarios) 


Catch variability was high for the highest value of op 
(scenarios 3, 20-21, and 27-29) and very high when op 
and pr were at their maxima (scenarios 35-37 and 52-53 
in Figure 6). In contrast, catch variability was insensi- 
tive to the other factors considered and to the choice of 
management strategy (Fig. 6). The probabilities of catch 
being less than mean catch and CPUE being less than 
mean CPUE were very high when op and pz were at their 
maxima (scenarios 35-37 and 52-53) (see Figure 7 for 
the linear choice and Supplementary Figure 6 [online only] 
for the nonlinear choice). The probability of not achieving 
mean CPUE was higher for HR30 (Fig. 7H, Suppl. Fig. 6H 


Annual proportional change in catch 


Annual proportional change in catch 


0.000 (1) 
1508 (4) 
0.049 (1) 
0.468 (4) 
47.4 (1) 
0.000 (1) 
16,315 (1) 
0.002 (1) 


0.000 (1) 
1732 (3) 
0.052 (3) 
0.140 (2) 
36.7 (2) 
0.153 (2) 
24,202 (2) 
0.075 (2) 


0.000 (1) 
1735 (2) 
0.055 (4) 
0.141 (3) 
36.4 (3) 
0.174 (3) 
24,444 (3) 
0.078 (3) 


0.000 (1) 
1770 (1) 
0.052 (2) 
0.101 (1) 
33.0 (4) 
0.408 (4) 
27,506 (4) 
0.196 (4) 


0.000 (1) 
1354 (4) 
0.048 (1) 
0.807 (4) 
139.6 (1) 
0.000 (1) 
4974 (1) 
0.037 (1) 


0.000 (1) 
1510 (2) 
0.050 (3) 
0.477 (2) 
103.8 (2) 
0.000 (1) 
7460 (2) 
0.396 (2) 


0.000 (1) 
1512 (1) 
0.052 (4) 
0.475 (1) 
103.5 (3) 
0.000 (1) 
7492 (3) 
0.398 (3) 


0.000 (1) 
1508 (3) 
0.049 (2) 
0.503 (3) 
96.7 (4) 
0.000 (1) 
7997 (4) 
0.535 (4) 


[online only]) than for HR10, HR15, and HR15U for all 
scenarios, with probability levels lower for the nonlinear 
choice than for the linear choice. 

Figure 8 and Supplementary Figure 7 (online only) show 
the probabilities of MMA being less than MMA,,,, and 
of MMB being less than MMB; for the linear and non- 
linear choices, respectively, when MMB.o;, was greater 
than MMB,, for the 5 HCRs. Spawning biomass and 
abundance declined even under HRO (i.e., F=0) for sce- 
narios 35-37 and 52-53 because of high recruitment 
variability for the base level of the steepness parame- 
ter (Fig. 8, A and F; Suppl. Fig. 7, A and F [online only]). 
The probabilities of MMB being less than MSST and of 
MMA being less than 0.25MMA,,,. were closer to zero for 
all scenarios and HCRs across simulations and for the 
entire projection period (results not shown). 

The median rebuilding time for all HCRs under the lin- 
ear choice increased for scenarios 35-37 and 52-53 
because of high recruitment variability (Fig. 9). As 
expected, HRO had the shortest rebuilding time. The 
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Figure 6 


Median catch variability in the directed pot fishery for golden king crab (Lithodes aequispinus) in the Aleutian Islands 
during the last 10 years of a 30-year projection period, which begins in 2018, for 53 scenarios of an operating model used to 
evaluate harvest control rules (HCRs). The stock is projected from 2 initial levels of abundance, measured in mature male 
biomass (MMB). Values are given in numbers without a unit for HCRs (A) HR10, (B) HR15, (C) HR15U, and (D) HR30, with 
the stock set initially at a healthy state (.e., MMB.0;:/MMB3,=1.55, where MMBoo1g is MMB in 2018 and MMB,, is 35% of 
the unfished level of MMB), and for HCRs (E) HR10, (F) HR15, (G) HR15U, and (H) HR380, with the stock set initially at an 
overfished state (i.e., MMB. ;./MMB,-=0.50). In the model used in this analysis, a linear relationship between catch per unit 
of effort and selected abundance is assumed. Data used in the model are for golden king crab in 1981-2018. For details about 


the HCRs, see Table 1. 


rebuilding time is short when values for steepness are 
higher (e.g., scenario 9) and is longer for HR30 than for 
HR10, HR15, and HR15U for most scenarios. The rebuild- 
ing times necessary to achieve MMA,,, are longer than 
those necessary to achieve MMB,, for all non-zero exploita- 
tion rates. The results for the nonlinear choice are similar 
(results not shown). 

The MMA and CPUE distributions for scenario 1 (with 
parameters estimated through the use of the assess- 
ment model) indicate that HR10, HR15, and HR15U per- 
formed better than HR30. Fewer instances of MMA being 
lower than MMA,,,. and more frequent values of higher 
CPUE were observed for HR10, HR15, and HR15U than 
for HR30 (Fig. 10). Catch distributions were similar for 
HR15 and HR380 and for the equilibrium state; size com- 
positions of crab in total catch as well as in retained catch 
did not differ significantly from those for HR15 and those 
for HR30 (Suppl. Fig. 8) (online only). This result occurred 
because the stock was projected from a healthy state 


(MMB>MMB,;;) and because the effective F under both 
HR15 and HR30 was dampened because of reduction in 
abundance over the years. 

The results for the conservation and economic perfor- 
mance criteria of HR15 and HR15U for the last 10 years of 
the 30-year projection period were very similar for the lin- 
ear and the nonlinear relationships of CPUE to abundance 
(Tables 4-6, Suppl. Tables 7—9 [online only]). When compared 
with HR10 and HR30, HR15 may optimize the balance 
between conservation and economic criteria (Tables 4-6). 
Conservation metrics were similar among all HCRs except 
for HR30 (nonlinear choice; Suppl. Table 7 [online only]). 
Probabilities of fishery closures were zero, and catch vari- 
ability was similar for all HCRs, yet relative stock status 
was lower under HR30. Relative to HR15, HR10 required 
substantially lower effort (~48% lower) to achieve the TAC, 
yet the TAC itself was lower (~15% lower); in contrast, 
HR30 required higher effort (~14% higher) to achieve only 
a marginally higher (~2% higher) TAC (Table 5). 
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Figure 7 
The probability of (A-D) catch and (E—H) catch per unit of effort (CPUE) being below their mean values for the period 
after crab stocks were rationalized, from the 2005-2006 fishing season through the 2018-2019 fishing season, in the 
directed pot fishery for golden king crab (Lithodes aequispinus) in the Aleutian Islands, by model scenario, under harvest 
control rules (HCRs) HR10, HR15, HR15U, and HR80. Estimates are based on the last 10 years of a 30-year projection 
period for 53 scenarios of an operating model used to evaluate HCRs. The stock is projected from an initial level of abun- 


dance, measured in mature male biomass (MMB): 1.55MMB,,, where MMB,, is 35% of the unfished level of MMB. In the 
model used in this analysis, a linear relationship between CPUE and selected abundance is assumed. For details about 


the HCRs, see Table 1. 


Discussion 


We have demonstrated the utility of the simplified MSE 
for evaluating HCRs through consideration of conserva- 
tion and economic trade-offs. Our results indicate that 
HR15 and HR15U are preferable to HR10 and HR30, 
given the desire to balance the trade-off between sus- 
tainability and economic viability, and this notion was 
supported under both the linear and nonlinear assump- 
tions about the relationship between CPUE and selected 
abundance. Although HR30 yielded the highest catch, it 
performed poorest in terms of other economic and con- 
servation criteria. Specifically, HR30 had the highest rel- 
ative probabilities of MMB being below MMB,, and of 
MMA being below the historical average MMA, had the 
lowest CPUE, and required substantially more effort to 
realize marginally higher (~2% higher) catches compared 
with HR10, HR15, and HR15U. Relative to HR10, HR30 
required 69% more effort to achieve 17% more catch 
(Table 5). Therefore, criteria beyond projected average 


TAC are important from an economic viewpoint because 
costs required to make excessive numbers of fishing 
trips when fishing effort is high may outweigh modest 
increases in TAC. 

Although we suggest that HR15 is the optimal HCR, 
given the trade-offs between conservation, catch, and 
catch stability, we acknowledge that the preferred HCR 
may differ depending on management and stakeholder 
priorities. Although HR10 yielded improved performance 
in terms of some conservation and economic criteria 
(e.g., higher CPUE, MMA, MMB, and reduced effort), 
it yielded lower catch compared with that from HR15, 
HR15U, and HR30. Our analysis is meant to provide 
managers and stakeholders with a tool to evaluate the 
trade-offs between various fishery management criteria 
relative to risk. 

The conservation and economic criteria were very simi- 
lar for HR15 and HR15U. Although the catch limit on 
legal-sized male abundance did not meaningfully affect 
the performance of the management strategy, it is likely 
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Figure 8 


The probability (A-E) of mature male abundance (MMA) being less than average MMA (MMA,,,.) and (F—J) of mature male 
biomass (MMB) being less than 35% of the unfished level of MMB (MMB,,;) for golden king crab (Lithodes aequispinus) in the 
Aleutian Islands during the period from the 1985-1986 fishing season through the 2018-2019 fishing season, by model scenario, 
under harvest control rules (HCRs) HRO, HR10, HR15, HR15U, and HR30. Estimates are based on the last 10 years of a 30-year 
projection period for 53 scenarios of an operating model used to evaluate HCRs. The stock is projected from an initial level of 
abundance: 1.55MMB,,;. In the model used in this analysis, a linear relationship between CPUE and selected abundance is 


assumed. For details about the HCRs, see Table 1. 


an important conservation component of the management 
strategy given the desire to ensure future recruitment. 
Realized exploitation rates on abundance of legal-sized 
male golden king crab is expected to be higher when pop- 
ulation abundance is on an increasing trend (i.e., when 
mature male recruits have yet to reach the legal size) 
because the exploitation rate is scaled to MMA. The maxi- 
mum exploitation rate on legal-sized male abundance pro- 
vides an additional level of protection against overfishing 
of legal-sized males in years when legal-sized male abun- 
dance is low relative to the abundance for the entire size 
range of mature males and is a commonly adopted step in 
the HCRs for other crab stocks in the Bering Sea and 
Aleutian Islands (e.g., red king crab in Bristol Bay; Pen- 
gilly and Schmidt’; Zheng et al., 1997). 


* Pengilly, D., and D. Schmidt. 1995. Harvest strategy for Kodiak 
and Bristol Bay red king crab and St. Matthew Island and 
Pribilof blue king crab. Alaska Dep. Fish Game, Spec. Publ. 7, 10 
p. [Available from website.] 


It is important to note that our simulations limited 
the catch in the directed fishery by the retained catch 
component of the ABC. As such, the more aggressive 
HCRs likely performed more conservatively (i.e., there 
was zero probability of exceeding OFL and ABC) than 
any of the HCRs that did not constrain the directed 
fishery catch below the retained catch component of the 
ABC. However, our simulations best approximate how 
management for crab stocks occurs in the North Pacific 
Ocean because TACs would not be set above estimated 
ABCs in practice. 

Incorporating uncertainty is a fundamental challenge 
in MSE. Parameters, such as natural mortality, catch- 
ability, growth, maturity, selectivity, and the stock— 
recruitment relationship are assumed to be correct and 
time invariant, and the projections ignore spatial and 
environmental variability (Somerton and Otto, 1986; 
Hollowed et al., 2001; Clark and Hare, 2002). However, 
the ability to use MSE to help achieve management goals 
depends on how well uncertainty in the system is repre- 
sented in simulations (Punt et al., 2016). Because it is 
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Figure 9 


Median time, in years, for a stock of golden king crab (Lithodes aequispinus) in the Aleutian Islands (A—E) to rebuild from an 
initial overfished state of 0.5MMB;;, where MMB; is 35% of the unfished level of mature male biomass, to a full MMB,, and 
(F—J) to rebuild from the corresponding mature male abundance to average mature male abundance under harvest control 
rules (HCRs) HRO, HR10, HR15, HR15U, and HR30, for the 53 scenarios of the operating model in which a linear relation- 
ship between catch per unit of effort and selected abundance is assumed. Data used in the model are for golden king crab in 


1981-2018. For details about the HCRs, see Table 1. 


computationally prohibitive to address all uncertainties 
by using a projection model, we considered a small set of 
scenarios focused on those uncertainties most likely to 
affect the performance of HCRs under the linear and 
nonlinear choices. 

Our findings indicate how projection results respond 
to changes in steepness, variation, and autocorrelation of 
the stock—-recruitment relationship, error in estimating 
MMB and MMA, and catch implementation error (Table 2, 
Suppl. Tables 1—5 [online only]). Although the values of the 
performance metrics differ between the linear and nonlin- 
ear choices, trends in HCR ranks were largely unchanged 
between these choices, indicating that our analysis for 
these choices is robust for evaluating policy trade-offs. The 
results presented here are based on approximate closed- 
loop simulations because, although errors in estimating 
MMB and MMA are considered, the full stock assessment 
is not simulated because of computational limitations, 
and our analysis is not a full MSE. Comparison of the full 
suite of scenarios (i.e., a range of contrasting parameter 
values) reveals the level of risk related to each source of 


uncertainty when relying on best estimates of parameters 
for decision-making. 

The projections in our study identify recruitment vari- 
ability as the most important factor determining the per- 
formances of the HCRs, yet understanding causes of 
recruitment fluctuations is a fundamental challenge in 
modeling crab population dynamics. The results of our 
simulations may underestimate recruitment variability 
or fail to capture the non-stationarity of the nature of 
recruitment, and such underestimation or failure may 
bias estimates of HCR performance. Well-defined stock— 
recruitment relationships are rare for crab and lobster 
species because the underlying physical and biological 
processes that influence larval survival to the juvenile 
stages are difficult to define (Wahle, 2003). For red king 
crab in Bristol Bay, recruitment trends are consistent 
with decadal climate shifts (Zheng and Kruse, 2003), 
indicating the importance of environmental factors 
(Zheng and Kruse, 2006). Nevertheless, because of uncer- 
tainties in the stock—recruitment relationship (or lack 
thereof) for golden king crab in the Aleutian Islands, 
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Figure 10 


Distributions of (A-E) mature male abundance (MMA), (F-I) retained catch, and (J—M) catch per unit of effort (CPUE), all 
given in frequency counts, of golden king crab (Lithodes aequispinus) in the Aleutian Islands under harvest control rules (HCRs) 
HRO, HR10, HR15, HR15U, and HR30. Distributions are from 1000 simulations for scenario 1 of the operating model in which a 
linear relationship between CPUE and selected abundance is assumed. The stock is projected from an initial level of abundance, 
measured in mature male biomass (MMB): a healthy state (i.e., MMB.o;</MMB,,=1.55, where MMB.o;, is MMB in 2018 and 
MMB,,; is 35% of the unfished level of MMB). The vertical thick dashed lines in the top row of panels indicate average mature 
male abundance (MMA,,,.), and the vertical thin dashed lines indicate 0.25MMA,,,. Data used in the model are for golden king 


crab in 1981-2018. For details about the HCRs, see Table 1. 


understanding the sensitivity of the performance of can- 
didate HCRs to changes in recruitment parameters will 
assist managers in decision-making. 

We feel that the approach of using a simplified MSE 
presented here is a fair balance between a robust analy- 
sis with reduced computational demands and a full MSE 
that can be applied to any stock that is hard to age and 
for which there are several candidate HCRs. Although 
we highlight the importance of recruitment variability, 
managers of other stocks should consider their man- 
agement goals, dynamics of the stock that they manage, 


uncertainties, and candidate catch policies. The ability 
to objectively evaluate conservation and economic trade- 
offs leads to a level of transparency between managers 
and fishermen that allows productive dialogue regarding 
emerging HCRs. 
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Table 6 


Decision matrix based on average policy ranks within each 
metric used to evaluate harvest control rules (HRO, HR10, 
HR15, HR15U, and HR30) for golden king crab (Lithodes 
aequispinus) in the Aleutian Islands, when a linear rela- 
tionship between catch per unit of effort and selected 
abundance is assumed in the operating model. Values are 
based on the last 10 years of the 30-year projection period, 
which begins with 2018. The stock is projected from 2 ini- 
tial levels of abundance, measured in mature male bio- 
mass (MMB): a healthy state (i.e., MMBoo)./MMB3,=1.55, 
where MMBoo13 is MMB in 2018 and MMB,, is 35% of 
the unfished level of MMB) and an overfished state (i.e., 
MMB 499/MMB,;=0.50). Values are average ranks within 
each metric, with ranks of the average metric ranks given 
in parentheses. The ranks of the catch metric correspond 
only to the long-term averages of retained catch; therefore, 
no average ranks were computed. For details about the 
harvest control rules, see Table 1. 


Catch 
Metric Conservation Catch stability 
Healthy 
HRO 1.00 (1) 
HR10 1.00 (1) 
HR15 1.40 (3) 
HR15U 1.40 (3) 
HR30 1.80 (5) 
Overfished 
HRO 1.00 (1) 
HR10 1.00 (1) 
HR15 1.50 (3) 
HR15U 1.75 (4) 
HR30 2.00 (5) 


1.43 (1) 
2.00 (2) 
2.86 (3) 
2.86 (3) 


1.43 (1) 
1.86 (2) 
2.29 (3) 
2.71 (4) 


suggestions on the projection simulations and stakehold- 
ers of the fishery for golden king crab in the Aleutian 
Islands for suggestions on management strategy scenar- 
ios. This manuscript is contribution PP—285 of the Com- 
mercial Fisheries Division of the Alaska Department of 
Fish and Game, Juneau, Alaska. 
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Abstract—Improved understanding of 
the seasonal distribution, habitat use, 
and fishery interactions of the common 
thresher shark (Alopias vulpinus) in the 
western North Atlantic Ocean (WNA) is 
required for future management. We 
compiled and analyzed 3478 fishery- 
dependent capture records in the WNA 
between 1964 and 2019 to examine 
dynamics by sex and life stage (i.e., 
young of the year, juvenile, and adult). 
Sharks were captured over a broad geo- 
graphic range from the Gulf of Mexico to 
the Grand Banks, primarily in continen- 
tal shelf waters shallower than 200 m. 
Seasonal north-south movements along 
the east coasts of the United States 
and Canada were observed for all life 
stages and both sexes, with individu- 
als generally occurring at more north- 
erly latitudes in the summer and more 
southerly latitudes in the winter. Dis- 
tinct areas of more frequent capture 
in fisheries were identified for all life 
stages throughout their range. Com- 
mon thresher sharks were observed in 
waters with sea-surface temperatures 
of 4-31°C, most commonly of 12—18°C. 
The results of this study will help to 
identify essential fish habitat for each 
life stage of common thresher sharks 
along the U.S east coast and to develop 
management measures for the WNA 
population. 
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The common thresher shark (Alopias 
vulpinus) is one of 2 species in the fam- 
ily Alopiidae that occurs in temperate 
and subtropical regions of the North 
Atlantic Ocean (Compagno, 2001). 
In the western North Atlantic Ocean 
(WNA), common thresher sharks are 
thought to represent a discrete popu- 
lation (Castro, 2011; Rigby et al., 2019) 
and occur from Cuba to Newfoundland, 
Canada (Compagno, 2001). There is 
currently no large-scale, directed com- 
mercial fishery for this species in the 
WNA, but common thresher sharks 
are often caught incidentally in pelagic 
longline and demersal gill-net fisheries 
and occasionally are retained for sale in 
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the United States (NMFS, 2019a). The 
common thresher shark is also the tar- 
get of an extensive recreational shark 
fishery from Virginia to Maine with 
the majority of sharks landed (NMFS, 
2019a) because their flesh is considered 
excellent (Compagno, 2001). 

There is limited and conflicting 
information on the status of the pop- 
ulation of common thresher sharks in 
the WNA. Results of studies examining 
pelagic longline logbook data collected 
since the late 1980s for thresher sharks 
(i.e., the common thresher shark and 
the bigeye thresher, A. superciliosus, 
combined) indicate that their relative 
abundance has declined by 63-80% 
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over the time series (Baum et al., 2003; Cortés et al., 2007). 
In contrast, pelagic longline fisheries observer data indi- 
cate that relative abundance had either stabilized at low 
levels (Baum and Blanchard, 2010) or increased (Cortés 
et al., 2007) during the 1990s into the early 2000s. Results 
of recent analyses of pelagic longline logbook and observer 
data for common thresher sharks from the period 1992— 
2013 also indicate a stabilized trend in relative abundance 
at low levels (Young et al., 2016) to an overall increasing 
trend when temperature profiles are incorporated into the 
analyses (Lynch et al., 2018). Relative abundance indices 
derived from data from 5 long-term recreational sportfish- 
ing tournaments based out of New York and New Jersey 
also indicate an overall increasing trend but with high 
annual variability (Young et al., 2016). 

Although recent analyses provide an optimistic out- 
look, the true extent to which these trends are represen- 
tative of the population of common thresher sharks in the 
WNA remains unknown. For example, although observer 
data are often considered more reliable than fishery log- 
book data, especially for bycatch species, the sample size 
for observer data is much smaller and interactions of 
individuals with the pelagic longline fishery are likely 
limited by the coastal and temperate distribution of this 
species (Young et al., 2016). Nonetheless, this species 
is currently listed as vulnerable by the International 
Union for Conservation of Nature (IUCN) in the IUCN 
Red List of Threatened Species, with a decreasing popu- 
lation trend assumed throughout its global range (Rigby 
et al., 2019). There are also no formal stock assessments, 
international management measures, or species-specific 
quotas for common thresher sharks in the WNA, and 
there are only limited biological and ecological data for 
the development of domestic and international manage- 
ment policies. 

Detailed information on seasonal distribution, habitat 
use, migration patterns, and population structure of com- 
mon thresher sharks in the WNA is limited. Although 
records from fishery-dependent catch data support a 
north-south movement pattern along the U.S. east coast 
(e.g., Castro, 2011; Natanson and Gervelis, 2013), pub- 
lished movement data are limited to fishery-dependent 
data from the recapture of 2 individuals from the 230 
common thresher sharks tagged by the National Marine 
Fisheries Service (NMFS) Cooperative Shark Tagging 
Program from 1962 through 2018 (tag return rate: 
~0.9%). Both recaptured sharks had a north-south move- 
ment pattern along the northeastern coast of the United 
States (Kohler and Turner, 2019). However, the records 
of these recaptures provide minimal insight into the spa- 
tial and temporal extent of the annual movements of this 
species, given the extended times at liberty (1533 and 
2934 d) and the fact that both recapture events occurred 
in June. Migratory routes are also poorly known (Castro, 
2011), and only general descriptions of distribution and 
habitat use by sex or life stage have been presented from 
historical catch data. 

Given the uncertainty over the status of the population 
of common thresher sharks in the WNA, there is a need 


to improve our understanding of this species’ life history 
and general susceptibility to capture in commercial and 
recreational fisheries. Improved knowledge is needed con- 
sidering that their life history parameters (i.e., relatively 
slow growth, late age to maturity, and low fecundity; 
Gervelis and Natanson, 2013; Natanson and Gervelis, 
2013) increase their susceptibility to population decline 
(Smith et al., 1998) and that ecological data on their tem- 
perature and habitat preferences are few. Accordingly, 
we compiled and summarized fishery-dependent data on 
catch of common thresher sharks throughout the WNA 
1) to document their seasonal distribution and habitat 
use by sex and life stage, 2) to better understand fishery 
interactions, and 3) to provide information that facili- 
tates the identification of essential fish habitat (EFH) by 
life stage in U.S. waters. 


Materials and methods 


At-sea observer and other catch data on common thresher 
sharks originated mainly from the NMFS, specifically 
from the Northeast Fisheries Observer Program, Pelagic 
Longline Observer Program of the Southeast Fisheries 
Science Center, Shark Bottom Longline Observer Pro- 
gram, Cooperative Shark Tagging Program of the Apex 
Predators Program, Cooperative Tagging Center of the 
Southeast Fisheries Science Center, and Large Pelag- 
ics Survey. Additional records were obtained from the 
Industry Surveys Database of Fisheries and Oceans 
Canada (H. Bowlby, personal commun.), which is used 
to archive at-sea observer information from commer- 
cial fisheries. Data from each source were subjected to 
editing and quality control a priori by the provider and 
represent all records of common thresher sharks avail- 
able in their respective databases through 2019. Data 
requested from each source include the date of capture, 
location of capture (i.e., latitude and longitude), shark 
sex, shark length (i.e., fork length [FL] and total length 
[TL]) and nature (i.e., measured or estimated) of the 
length, shark weight (i.e., whole [or round] and dressed) 
and nature (i.e., measured or estimated) of the weight, 
gear type used for capture, and the observed sea-surface 
temperature (SST, in degrees Celsius) at the time and 
location of capture. 

Prior to analysis, records with no geographical informa- 
tion were removed, and a series of steps were taken to stan- 
dardize and validate length and weight data among the 
various sources. All capture records for which a measured 
FL was provided were considered reliable and retained. 
Records with a measured or estimated TL or whole weight 
were converted to an estimated FL by using the appropri- 
ate conversion equation (TL to FL: coefficient of multiple 
determination [R?]=0.83; whole weight [in kilograms] to 
FL: R?=0.93; NMFS’). Reported weights were also used to 


! NMFS (National Marine Fisheries Service). 2020. Unpubl. data. 
Apex Predators Program, Northeast Fish. Sci. Cent., Natl. Mar. 
Fish. Serv., NOAA, 28 Tarzwell Dr., Narragansett, RI 02882. 
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reconcile any lengths that appeared to have been inaccu- 
rately labeled (e.g., labeled as TL when the weight was 
more consistent with a fish near the same FL). Length 
data for any record that had a measured or estimated FL 
below the published size of near-term embryos (50 cm FL; 
Gervelis and Natanson, 2013) or markedly above the esti- 
mated maximum length (350 cm FL; Cailliet et al., 1983) 
were considered unreliable and excluded from ontoge- 
netic analyses. Lastly, all common thresher sharks with 
a measured or estimated FL were assigned to 1 of 3 life 
stages, on the basis of published estimates of length at age 
(Gervelis and Natanson, 2013) and length at 50% matu- 
rity (Natanson and Gervelis, 2013) for individuals from 
the WNA: young of the year (YOY), if <100 cm FL; juve- 
nile, for males 101—187 cm FL and for females 101—215 cm 
FL; and adult, for males >188 cm FL and for females and 
sharks of unknown sex >216 cm FL. 


Distribution analysis 


Spatial and temporal patterns in presence and fish- 
ery interactions of common thresher sharks by sex, life 
stage, and season were qualitatively analyzed following 
the general approach of Curtis et al. (2014). All analyses 
were performed in R (vers. 4.0.3; R Core Team, 2020). 
Maps and figures were created by using the tidyverse 
collection of packages (vers. 1.3.0; Wickham et al., 
2019) and the sf package (vers. 0.9.6; Pebesma, 2018) in 
R. Seasons were classified as winter (January—March), 
spring (April-June), summer (July-September), or fall 
(October-December). 

Because of the limitations of using presence-only data 
when detectability is unknown, location data over the full 
time series (1964-2019) were plotted by season, both in 
aggregate and by sex and life stage, to assess distribution 
patterns. To meet confidentiality requirements for com- 
mercially derived data, latitudes and longitudes (under 


the World Geodetic System 1984) were aggregated in a 
raster grid (0.5° x 0.5°) spanning the WNA, including the 
Gulf of Mexico. Grid cell counts (number of captures) were 
natural log transformed for plotting to increase contrast 
among small values. 


Habitat use 


Fishery-dependent records with a recorded SST were 
aggregated by month and summarized by using box 
plots to examine monthly and seasonal trends. To visu- 
ally examine distribution of common thresher sharks in 
relation to typical seasonal SST conditions in the WNA, 
grid cells were plotted over NMFS Southwest Fisheries 
Science Center monthly composite climatologies (data 
set ID: erdAGsstamday_LonPM180; available from the 
ERDDAP server at website) averaged over the most 
recent years for which data were available (2009-2016); 
62% of all capture events occurred over this period. The 
SST data were downloaded by using the rerddap package 
(vers. 0.6.5; Chamberlain, 2019) in R. To explore any asso- 
ciation with bathymetry, the local depth was assigned to 
each record by using data from the NOAA ETOPO1 1 
Arc-Minute Ocean Relief Model (NOAA National Geo- 
physical Data Center, model available from website, 
accessed June 2019; Amante and Eakins, 2009) in the R 
package marmap (vers. 1.0.5; Pante and Simon-Bouhet, 
2013). Histograms were used to examine SST and depth 
occupancy by life stage. 


Results 


A total of 3478 fishery-dependent records of the capture of 
common thresher sharks were compiled between 1964 and 
2019, representing 1035 males, 1039 females, and 1404 
sharks of unknown sex (Tables 1 and 2). The number of 


Table 1 


Source, number, and time period for fishery-dependent capture records com- 
piled to examine spatial distribution and habitat use of common thresher 
sharks (Alopias vulpinus) in the western North Atlantic Ocean between 1964 
and 2019. DFO=Department of Fisheries and Oceans; NMFS=National Marine 
Fisheries Service; SEFSC=Southeast Fisheries Science Center. 


Data source 


DFO Canada: Industry Surveys Database 
NMES Cooperative Shark Tagging Program 


NMFS Large Pelagics Survey 


NMES Shark Bottom Longline Observer Program 
NOAA Northeast Fisheries Observer Program 
NOAA Southeast Pelagic Observer Program 
NOAA SEFSC Cooperative Tagging Program 


Total 


Number of Time 
records period 


299 1982-2019 
1964-2019 
2002-2018 
1996-2019 
1990-2019 
1992-2017 
1978-1996 
1964-2019 
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Table 2 


Summary of fishery-dependent data by life stage and sex as well as the range of associated sea-surface 


temperatures (SST) and depths from capture records 


for common thresher sharks (Alopias vulpinus) 


in the western North Atlantic Ocean between 1964 and 2019. Ranges and means, with standard devi- 
ations in parentheses, are presented for SST and depth. YOY=young of the year; unknown=unknown 


life stage or sex. 


Number of records 


Life stage Male Female Unknown 


YOY 425 432 279 
Juvenile 338 376 369 
Adult 182 83 158 
Unknown 90 148 598 
Overall 1035 1039 1404 


Bottom longline Haul seine 
Gear type 


Percentage of records 


TOKTRODOTANTOORDOOTAMTOORD 
OOOMKEKEN DODDODDDDOHOHMWWOOOOMAO 
foptopToploploploploploplopioploptoploploploptopioplopioploptopioplopiopieo) 


AAS 


Year 


Figure 


SST (°C) Depth (m) 


Range Mean Range Mean 


6-26 16 (3) 1—2077 21 (86) 
6-31 17 (4) 1—4535 262 (556) 
5-30 18 (4) 1-5401 414 (902) 
4-30 18 (4) 1-5427 1197 (1581) 
4-31 17 (4) 1-5427 434 (1039) 


Otter trawl | Purse seine 


Gill net Mid-water trawl | Pelagic longline | Rod and reel 


syJeus Jo Jequinn 


DONO OOOO NONIOOR OGY 
HOOOCOCOOCOCOO— 

DOOSOOSCOSOSOOOOOOOOOO 
TANNNAANNANNANANANNNNANAOAOAN 


1 


Percentage of the total number of capture records from fishery-dependent data, by gear type, and 
the total number of records, by year, for common thresher sharks (Alopias vulpinus) in the western 
North Atlantic Ocean from 1964 through 2019. The shade or pattern of each bar indicates gear 


type, and the solid line indicates number of records. 


records obtained per year ranged from 1 to 538 (mean: 76 
[standard deviation (SD) 109]) with 48.7% of all records 
occurring between 2015 and 2019 (Fig. 1). Lengths were 
available for 3072 (88.3%) individuals, including 1720 
(49.4%) sharks that were measured and 1352 (38.9%) 
sharks whose lengths were estimated either at the time of 
capture or by a length conversion. Measured FLs ranged 
from 59 to 320 cm (mean: 125 cm [SD 48]), and estimated 


FLs ranged from 60 to 330 cm (mean: 160 cm [SD 56]) 
(Fig. 2). Five length records were eliminated because of 
biologically implausible length estimates or measure- 
ments. Life stage was assigned to 2642 (76.0%) sharks of 
which 1136 (32.7%) were YOY, 1083 (31.1%) were juvenile, 
and 423 (16.0%) were adult (Table 2, Fig. 2). Sex was noted 
for only 21 of the sharks captured in Canada; therefore, 
any comparisons among sexes and juvenile and adult life 
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Figure 2 


Length-frequency histogram for 3072 common thresher sharks (Alopias vulpinus), representing 
males, females, and those with unknown sex for which a measured or estimated length was avail- 
able, captured from 1964 through 2019 in the western North Atlantic Ocean. The classification of 
life stages, based on sizes at 50% maturity from Natanson and Gervelis (2013), is provided at the 
top of the figure for males and females. YOY=young of the year. 


stages almost exclusively represent sharks captured in 
US. waters. 


Gear interactions 


Common thresher sharks were captured with a wide 
range of gear types, including commercial bottom long- 
line, pelagic longline, gill net (sink and floating), haul 
(beach) seine, purse seine, (bottom) otter trawl, mid-water 
trawl, and recreational rod and reel (Fig. 1). Of these 
gears, gill nets (50.0%; primarily sink gill nets) and 
pelagic longlines (27.1%) accounted for the majority of 
capture records, followed by rod and reel (10.4%) and otter 
trawl (9.2%). All other gears accounted for no more than 
1.4% of the total number of records. The most common life 
stage of sharks captured by gill nets was YOY (88.4%), but 
YOY sharks were also caught by otter trawls, haul seines, 
and rod and reel (Fig. 3). Juveniles were captured by all 
gear types but were most commonly taken in gill nets 
(43.1%) and pelagic longlines (33.3%). Adults were cap- 
tured by every gear type except haul seine, but they most 
commonly were taken by pelagic longline (39.2%) and rod 
and reel (26.7%). 


Seasonal distribution 


Common thresher sharks were captured over a broad 
geographic range extending from 23.9°N to 46.3°N and 


from 42.8°W to 95.1°W (Fig. 4). The majority of records 
(number of sharks [n]=2751; 79.1%) were from conti- 
nental shelf waters shallower than 200 m; however, 
numerous individuals (n=727; 20.9%) were encoun- 
tered in deeper (>200 m) offshore waters of the WNA. 
Records were available from all months and seasons 
with the greatest number of them from the winter 
(Table 2). The increase in capture records during win- 
ter was partially due to a large number of YOY and 
juveniles being taken in 2018 by gill-net fisheries off 
North Carolina. 

Throughout the year, north-south changes in distribu- 
tion along the east coasts of the United States and Can- 
ada were noted for common thresher sharks over all life 
stages. In general, individuals occurred at more northerly 
latitudes in the summer and more southerly latitudes in 
the winter (Fig. 4). During the winter, sharks were primar- 
ily encountered off the coast of North Carolina, both in the 
vicinity of Cape Hatteras and in waters near the edge of 
the continental shelf. A considerable number of specimens 
were also caught in continental shelf waters off the east 
coast of Florida during the winter. Sharks were distrib- 
uted over the broadest area during the spring, with the 
center of distribution from North Carolina to Long Island, 
New York. During the summer, sharks were captured most 
commonly in continental shelf waters off the mid-Atlantic 
states, southern New England, and the Gulf of Maine, and 
sharks were rarely encountered south of Maryland. In the 
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Figure 3 


Number of common thresher sharks (Alopias vulpinus) 
caught in the western North Atlantic Ocean, based on 
compiled fishery-dependent data from 1964 through 2019. 
Sharks were captured with the following gear types: gill 
net (GN), pelagic longline (PLL), otter trawl (OT), rod 
and reel (RR), haul (beach) seine (HS), mid-water trawl 
(MWT), bottom longline (BLL), and purse seine (PS). Note 
the difference in scale between the y-axes. YOY=young of 
the year. 


fall, 2 centers of distribution were evident: one from North 
Carolina to New Jersey and another off the Scotian Shelf 
in Canada. 


Young of the year Common thresher sharks classified 
as YOY were distributed almost exclusively (99.1%) in 
continental shelf waters north of 33.5°N (Fig. 5). Only 6 
individuals were captured south of 33.5°N, including a 
single record from the Gulf of Mexico (27.8°N, 88.6°W), 
which may have been a misidentification. Sharks of this 
life stage migrated along the coast in a north-south 
pattern, with overwintering grounds off Cape Hatteras. 
During the summer, sharks were most commonly taken 
off New Jersey and Long Island, New York, but their 
distribution ranged as far north as Cape Cod, Massa- 
chusetts. During the spring and fall, their distribution 
spanned between these summer and overwintering 
grounds. 


Juveniles and adults Because of the inability to assign a 
life stage to all recorded common thresher sharks and 
the limited data from adults (Table 2), composite distri- 
bution maps were made for the juvenile and adult life 
stages (Fig. 6). The lack of sex data for sharks captured in 
Canada also precluded a thorough assessment of juvenile 


and adult distribution in waters of Canada, particularly 
along the Scotian Shelf. The winter distribution of these 
life stages occurred off the coast of North Carolina and 
the center of summer distribution extended from the 
mid-Atlantic states north to the Gulf of Maine. The distri- 
bution patterns of both juvenile and adult sharks are con- 
sistent with a seasonal, north-south distribution pattern 
along the east coasts of the United States and Canada, 
with extensive movements between the winter and sum- 
mer habitat areas. 


Males and females In general, the overall and seasonal 
distributions of male and female common thresher 
sharks were similar and followed the patterns previously 
described for this species and for each life stage (Fig. 7). 


Habitat use 


Sea-surface temperature at the time and location of cap- 
ture was recorded for 2764 common thresher sharks and 
ranged from 4°C to 31°C (mean: 17°C [SD 4]) (Table 2); 
78.2% of all records were for sharks captured in waters 
with temperatures of 10—22°C (Fig. 8). Young-of-the-year 
sharks were associated with the narrowest temperature 
range (6—26°C); however, each life stage and both sexes 
were generally associated with a similar range of SSTs 
(males: 5-30°C [mean: 17°C (SD 4)]; females: 6—30°C 
[mean: 17°C (SD 4)]). Monthly and seasonal fluctuations 
in SST at capture were evident for all life stages (Fig. 9). 
Sharks were taken in areas with depths ranging from 
1 to 5427 m (mean: 433 m [SD 1032]), with 43.7% of all 
records from waters 10-25 m deep (Table 2, Fig. 8). In 
particular, 98.4% of YOY were caught in depths from 1 
to 50 m, and 61.7% of juveniles were caught in depths 
from 10 to 100 m. Males and females were taken in 
areas with depths from 1 to 5427 m (mean: 165 m [SD 
533]) and from 1 to 5401 m (mean: 146 m [SD 497]), 
respectively. 


Discussion 


The results of our analysis of over 50 years of fishery- 
dependent data advance the scientific understanding of 
the distribution and habitat use of common thresher 
sharks in the WNA and will improve our ability to assess 
and sustainably manage the population of this species 
nationally and internationally. Although trends of relative 
abundance for common thresher sharks in the WNA 
appear to have stabilized or slightly increased in recent 
years, current levels are still well below historic values of 
abundance (Young et al., 2016; Lynch et al., 2018) and fish- 
ery interactions are increasing across some gear types 
(Fig. 1) (Gervelis and Natanson, 2013; NMFS, 2019a). New 
insights into habitat use by sex and life stage will also 
promote the identification of EFH in U.S. waters. Accurate 
EFH designations allow potential effects on species in 
areas of proposed development or other activities to be 
considered before approval and provide a basis for 
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Figure 4 


Distribution of common thresher sharks (Alopias vulpinus) in the western North Atlantic Ocean by season during 1964-2019 
and in relation to average sea-surface temperature (SST) climatologies for the period 2009-2016. The number of capture records 
(Rec) in each grid cell (0.5° x 0.5°) is natural log transformed. Records from the Gulf of Mexico are presented in the inset. The 


200-m depth contour delineates the boundary of the continental shelf. 


evaluating time or area closures to assist in the recovery 
of overfished populations. 

Analyses of presence-only data are a useful and 
cost-effective way to understand the distribution and 
habitat use of wide-ranging species for which limited 
species-specific survey data exist (Elith and Leathwick, 
2009; Curtis et al., 2014). Their main limitation is the 
potential for bias due to spatial and temporal variabil- 
ity in fishing effort or in catchability by both gear type 
and life stage (Pearce and Boyce, 2006; Curtis et al., 
2014). For example, common thresher sharks are the 
target of a relatively large, seasonal recreational shark 
fishery that operates from Virginia to Maine (Gervelis 
and Natanson, 2013), and >99% of the observed rod- 
and-reel capture events occurred over this geographic 
area. However, commercial fisheries data would be 
expected to be less biased, given that common thresher 
sharks are caught almost exclusively as bycatch in the 


WNA (Castro, 2011; Young et al., 2016) and that effort 
in trawl, gill-net, and pelagic longline fisheries is rela- 
tively consistent, widespread, and year-round along the 
U.S. east coast (Beerkircher et al., 2004; Guiet et al., 
2019). Because 94.9% of records originated from a vari- 
ety of commercial fisheries that operate throughout the 
range of this species, the overall and seasonal trends 
we report herein may largely reflect animal distribu- 
tion. Nonetheless, spatial and temporal concentration 
in fishing effort may have influenced the number of 
captured sharks by season; for example, the large num- 
ber of YOY observed off North Carolina during the win- 
ter and spring. 


Seasonal distribution and habitat use 


Our results are consistent with those of previous reports 
that indicate that the common thresher shark is primarily 
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Figure 5 


Distribution of young-of-the-year common thresher sharks (Alopias vulpinus) by season in the 
western North Atlantic Ocean between 1964 and 2019. The number of capture records (Rec) in 
each grid cell (0.5° x 0.5°) is natural log transformed. One record of a capture in the Gulf of Mexico 
is not presented because it may have been a misidentification. The 200-m depth contour delineates 


the boundary of the continental shelf. 


distributed at depths <200 m in the continental shelf 
waters in the WNA but that its range seasonally extends 
farther offshore, closer to the northern edge of the Gulf 
Stream along the Scotian Shelf and Grand Banks (Fig. 4) 
(Castro, 2011). Although higher rates of observer coverage 


in coastal fisheries (NMFS, 2019b) would contribute to 
more frequent capture of common thresher sharks in con- 
tinental shelf waters, their limited occurrence in the 
extensive and fairly well-monitored pelagic longline fish- 
eries in offshore waters of the United States and Canada 
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Figure 6 


Distribution of juvenile and adult common thresher sharks (Alopias vulpinus) by season in the western North Atlantic Ocean 
between 1964 and 2019. The number of capture records (Rec) in each grid cell (0.5° x 0.5°) is natural log transformed. The lack 
of records in the maps for waters of Canada is a result of missing data on sex of sharks and the subsequent inability to assign 
a life stage to sharks captured in Canada. Records from the Gulf of Mexico are presented in the inset. The 200-m depth contour 


delineates the boundary of the continental shelf. 


(Beerkircher et al., 2004; Hanke et al.”) indicates that this 
species does not occur as frequently in deep, offshore waters. 
Capture of common thresher sharks has been recorded spo- 
radically throughout tropical offshore waters, including in 
the Gulf of Mexico. However, species identification could not 
be verified, and it is possible that some capture records rep- 
resent misidentifications of bigeye thresher, the only other 
alopiid that occurs in the WNA (Compagno, 2001). Of note, 
misidentification would be most likely in the pelagic long- 
line records from offshore areas in the Gulf of Mexico, where 
the distribution of bigeye threshers is well-documented 


? Hanke, A. R., I. Andrushchenko, and G. Croft. 2012. Observer 
coverage of the Atlantic Canadian swordfish and other tuna 
longline fishery: an assessment of current practices and alter- 
native methods. DFO Can. Sci. Advis. Secr. Res. Doc. 2012/049, 
84 p. [Available from website.] 


(Fernandez-Carvalho et al., 2015) and the species tends to 
be encountered 4 times more frequently than the common 
thresher shark (Young et al., 2016). 

Trends in distribution are consistent with a seasonal, 
north-south movement pattern in the WNA, with distri- 
bution focused shoreward of the continental shelf edge in 
depths <200 m along the east coasts of the United States 
and Canada. This distribution and movement pattern 
is similar to those of many other highly migratory species 
that inhabit the temperate WNA, such as the shortfin 
mako ([surus oxyrinchus) (Vaudo et al., 2017) and bluefin 
tuna (Thunnus thynnus) (Galuardi et al., 2010), except 
with no evidence of frequent movements to the tropical 
Atlantic Ocean. It also aligns with the seasonal distribu- 
tion and movements of common thresher sharks in the 
eastern Pacific Ocean, where northward (spring and sum- 
mer) and southward (fall) movements occur between Baja 
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Figure 7 


Distribution of male and female common thresher sharks (Alopias vulpinus) and of both sexes combined in the western North 
Atlantic Ocean between 1964 and 2019, by season and aggregated in grid cells (0.5° x 0.5°). The lack of capture records in maps 
for waters of Canada is a result of missing data on sex of sharks. Records from the Gulf of Mexico are presented in the inset. The 
200-m depth contour delineates the boundary of the continental shelf. 


California, Mexico, to as far north as British Columbia, 
Canada (Hanan et al., 1993; Cartamil et al., 2016; Kinney 
et al., 2020). The distribution of common thresher sharks 
in the WNA was more contracted during the winter than 
in any other season, with individuals being most com- 
monly taken in relatively small areas off the coast from 
Virginia to North Carolina (outside of the barrier islands) 
and off the east coast of Florida. 

The seasonal distribution of common thresher sharks 
changes markedly over ontogeny. Parturition for this spe- 
cies occurs from May through August in the mid-Atlantic 
states and in southern New England (Natanson and 
Gervelis, 2013; Young et al., 2016), consistent with the ele- 
vated occurrence of YOY throughout these areas during 
the summer months (Fig. 5). In the winter, YOY were taken 
almost exclusively in depths <50 m off the coast of North 
Carolina from Oregon Inlet to Cape Lookout. These results 
align with previous hypotheses that inshore waters from 


North Carolina to Massachusetts support nursery habitat 
for YOY and juvenile common thresher sharks (Natanson 
and Gervelis, 2013) and indicate that specific areas off 
New York, New Jersey, and Cape Hatteras (Fig. 5) may 
be of particular importance to individuals in these life 
stages. Investigation into why YOY occur more frequently 
in these areas than in others was beyond the scope of this 
study; however, their presence may be related to environ- 
mental preferences (see subsequent paragraphs) or prey 
distribution. 

The distribution of common thresher sharks expands 
and individuals move into deeper waters of the WNA as 
they grow to maturity. A similar range expansion has been 
documented for juveniles and adults of this species off the 
U.S. west coast (Smith et al., 2008; Kinney et al., 2020) 
and may be related to changes in diet or physiological lim- 
itations imposed by environmental variation (Knip et al., 
2010). Note that the inability to assign a life stage to 
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female sharks in late summer (senior 
author and D. Bernal, unpubl. data). 
Sea-surface temperatures recorded 
at the time of capture indicate that the 
common thresher shark occurs over a 
wide range of temperatures but is most 
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Figure 8 


Percentage of the total number of capture records from fishery-dependent data 
of common thresher sharks (Alopias vulpinus) in the western North Atlantic 
Ocean between 1964 and 2019 by (A) sea-surface temperature (SST) and 
(B) depth for each life stage: young of the year (YOY), juvenile, adult, and 


unknown. 


24.0% of the captured sharks and the low number of 
adults (i.e., 423 sharks) hindered our ability to describe 
the distribution and habitat use of sharks in these life 
stages, particularly in Canada. However, the clear associ- 
ation of YOY with shallow waters (depths <50 m) indi- 
cates that records from deep, offshore waters of Canada 
were of juveniles and adults. Difficulties measuring or 
estimating the total length of common thresher sharks 
and variability in size at maturity (Natanson and Gervelis, 
2013) may have also led to the misclassification of life 
stage for some individuals and confounded our descrip- 
tions of distribution for each life stage. 

Contrary to research in other regions (Moreno et al., 
1989; Smith et al., 2008; Kinney et al., 2020), our results 
indicate that there is no strong evidence of sex-specific spa- 
tial or depth segregation of the population in the WNA. Seg- 
regation by sex tends to be most apparent in the adult life 
stage (Speed et al., 2010), and it is possible that we did not 
have sufficient information on this life stage for trends to be 
apparent. There is some evidence of sexual segregation of 
common thresher sharks in southern New England based 
on observations of large aggregations of predominantly 


commonly associated with temperatures 
from 10°C to 22°C. Previously, this spe- 
cies has been observed in SSTs ranging 
from 16.5°C to 19.8°C off the east coast of 
Florida (Castro, 2011), from 8°C to 28°C 
near the Republic of the Marshall Islands 
(Cao et al., 2011), and from 9°C to 21°C 
off the U.S. west coast (Cartamil et al., 
2016). Preferred temperatures have been 
reported as 18—20°C (Cao et al., 2011) and 
14-17°C (Cartamil et al., 2016), which 
are relatively similar to the temperature 
range over which the greatest number of 
capture events were recorded during this 
study (i.e., 12-18°C; Fig. 8). The overall 
range and mean SST at capture was also 
similar for all life stages (Table 2); how- 
ever, YOY were observed over the narrow- 
est temperature range. This result may 
be due to the lower capacity for metabolic 
heat retention in small individuals, which 
have a higher surface-to-volume ratio 
and therefore a decreased capacity to 
retain heat through regional endothermy 
(Bernal and Sepulveda, 2005). 

The monthly and seasonal patterns in 
SST at capture indicate that tempera- 
ture may be one of the key determinants 
of distribution and migration of common 
thresher sharks in the WNA. In general, 
captured sharks rarely occurred north of 
~37°N in the winter and south of ~37°N in the summer (Fig. 
4), the seasons during which the lowest and highest mean 
SSTs were observed, respectively (Fig. 9). This pattern is 
typical of migratory species that inhabit coastal waters of 
the WNA. By contrast, temperature was not one of the main 
factors driving movements and distribution off the U.S. 
west coast, with more broad-scale extrinsic seasonal factors 
(e.g., North Pacific Gyre Oscillation) being more influential 
(Kinney et al., 2020). Of note, the trends reported herein 
solely represent SST at the time and location of capture and 
do not necessary represent the temperatures that sharks 
experience at depth, which are likely the main drivers of 
movement and distribution. This is supported by Lynch 
et al. (2018), who demonstrated that water temperature 
profiles influenced catch rates for common thresher sharks 
in the pelagic longline fishery that operates off the Atlantic 
coast of the United States. 


Fishery interactions and management implications 


Capture records were available for common thresher 
sharks of all life stages; however, YOY was the most 
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Figure 9 


Box plot of sea-surface temperatures (SST) recorded, by month and season, 
as part of capture records from fishery-dependent data of common thresher 
sharks (Alopias vulpinus) in the western North Atlantic Ocean between 1964 
and 2019. The numerals above the boxes indicate the number of records for 
each month. The upper and lower parts of each box represent the first and 
third quartiles (the 25th and 75th percentiles), and the horizontal white line 
is the median. The whiskers that extend above and below the box correspond 
to 1.5 times the interquartile range, and the black circles represent values 


outside this range. 


commonly observed life stage, possibly because the 
restricted distribution of YOY in depths <50 m. There are 
no minimum commercial size limits for this species in 
the U.S. waters of the Atlantic Ocean and a review of the 
data provided by the Northeast Fisheries Observer Pro- 
gram indicate that many YOY are retained for sale. For 
example, between January and April 2018, at least 367 
YOY were captured by gill nets during observed trips in 
the vicinity of Cape Hatteras. Assuming an average 
brood size of 3.7 young (Natanson and Gervelis, 2013), 
these captured sharks represent the reproductive output 
of approximately 100 adult females. Given the bi- or tri- 
ennial reproductive cycle of this species in the WNA 
(Natanson and Gervelis, 2013) and the fact that the total 
number of YOY removed in the fishery may be many 
times greater than what was observed, such removals 
may represent a significant loss in recruitment in the 
WNA. However, because of their small size, landed 
weight of YOY common thresher sharks may not consti- 
tute a large percentage of the total quota for pelagic 
sharks in the Atlantic Ocean (i.e., 488.0 metric tons in 
dressed weight; NMFS, 2006). 

Our results provide new insight into the environmental 
conditions and habitats used by common thresher sharks 
for each life stage and will assist managers in the designa- 
tion of EFH by life stage off the U.S. east coast, as required 


Nov Dec 


under the Magnuson-Stevens Fishery 
Conservation and Management Act. 
Because of insufficient data on individ- 
ual life stages, all life stages currently 
share a single EFH that spans conti- 
nental shelf waters from Cape Lookout 
to Georges Bank and from Cape Ann, 
Massachusetts, to Bar Harbor, Maine 
(NMFS, 2017). However, our results 
provide information that should be suf- 
ficient to differentiate EFH for YOY. 
Although the extent of fishery-dependent 
data compiled for YOY common thresher 
sharks ranges from the Gulf of Maine 
to off South Carolina (consistent with 
Natanson and Gervelis, 2013), our anal- 
yses clearly identify areas of high levels 
of interaction off New York and New 
Jersey during the summer and off Cape 
Hatteras during the winter. Additional, 
direct information on the residency and 
movements of YOY common thresher 
sharks in these areas, however, is needed 
to confirm whether they warrant desig- 
nation as nursery areas (Heupel et al., 
2007). Additional fishery-independent 
(e.g., tagging and survey) data are also 
needed to better define juvenile and 
adult EFH and to confirm the broader 
conclusions regarding distribution and 
habitat use of common thresher sharks 
presented herein. 
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All submissions are subject to a double-blind review 
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Plagiarism and double publication are considered 
serious breaches of publication ethics. To verify the orig- 
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plagiarism-detection software. 

Manuscripts must be written in English; authors 
whose native language is not English are strongly advised 
to have their manuscripts checked by English-speaking 
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addresses, and the senior author’s email address. 
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not review the methods of the study or list the contents of 
the paper. Because abstracts are circulated by abstract- 
ing agencies, it is important that they represent the 
research clearly and concisely. 


General text must be typed in 12-point Times New 
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Stone et al. (2011).” 
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following: 


° Aging 
For our journal, the word aging is used to mean both 
age determination and the aging process (senescence). 
Authors should make clear which meaning is intended 
where ambiguity may arise. 


e Fish and fishes 
The plural of the word fish (a collective noun that 
implies individuals without regard to species) is fish. 
Example: The fish were collected by trawl net. 
Example: The numbers of fish collected that season 
were less than the numbers from previous years. 


The plural for fish species is fishes (a contrived plural 
used by taxonomists to mean several or more fish spe- 
cles) or one can use fish species (which is preferred in 
this journal for clarity across disciplines). 
Example: The fishes of Puget Sound [biodiversity is 
implied] or 
Example: The fish species of Puget Sound [preferred 
plural for clarity across disciplines]. 


e Crab and crabs, squid and squids, ete. 
The plural of the word crab (i.e., many individuals 
without regard to species) is crab. 
Example: The crab were sorted by weight. 
Example: Many red king crab were dying [Many 
individuals of one species of crab.] 


The plural of crab species is crabs (a word used by tax- 
onomists) or crab species (the latter is preferred in this 
journal for clarity). 

Example: These crabs were selected for treatment. 
[Different crab species are implied.] 
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meaning. ] 

Example: Snow crabs are found throughout the 
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Example: Two species of snow crab are found through- 
out the North Pacific Ocean and Bering Sea. [Pre- 
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ment. [Preferred word choice for clarity.] 
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We use fisherman and fishermen in this journal not fisher 
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e The definite article with common names of species 
When the singular common name of a species rep- 
resents the entire class or group to which it belongs, 
use the definite article. 

Example: Only one species of the genus Salmo is 
found in the Atlantic Ocean—the Atlantic salmon 
(Salmo salar). 

Example: The sonic emissions of the bottlenose dol- 
phin are complex. 


For plural common names, this rule does not apply. 
Example: Chinook salmon are found throughout the 
Pacific Ocean. 
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ate and tropical waters. 


e Sex 
For the meaning of male and female, use the word sex, 
not gender. Do not write “fish were sexed.” Write, “sex 
was determined.” 


e Participles 
As adjectives, participles must modify a specific noun 
or pronoun. 

Example: Using mark-recapture methods, these sci- 
entists determined the size of the population. [Correct. 
The participle wsing modifies the word scientists.] 

Example: These scientists, based on the collected 
data, concluded that the mortality rate of these fish 
had increased. [Incorrect. The scientists were not 
based on the collected data.] 

Example: These scientists concluded, on the basis of 
collected data, that the mortality rate of these fish 
had increased. [Correct. The offending participle has 
been eliminated and an adverbial phrase modifies 
the verb concluded.] 
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Equations and mathematical symbols should be set 
from a standard mathematical program (MathType or 
Equation Editor). Equations formatted in LaTex are 
not acceptable. For mathematical symbols in the gen- 
eral text (a, x2, ™, +, etc.), use the symbols provided by 
the MS Word program and italicize all variables, except 
those variables represented by Greek letters and the 
superscript and subscript parts of variables and expres- 
sions. Do not use photo mode when creating these 
symbols in the general text, and do not cut and paste 
equations, letters, or symbols from a different software 
program. 

Number equations (if there are more than one) for 
future reference by scientists; place the number within 
parentheses at the end of the first line of the equation. 


Literature cited section comprises published works and 
those accepted for publication (in press) in peer-reviewed 
journals. Follow the name and year system for citation for- 
mat in this section (i.e., citations should be listed alpha- 
betically by the authors’ last names, and then by year if 
there is more than one citation by the same author. A list 
of abbreviations for citing journal titles can be found on 
our website. 

Authors are responsible for the accuracy and com- 
pleteness of all citations. Avoid the use of multiple 
citations when a single citation sufficiently supports a 
statement; cite the work that first reported the informa- 
tion that supports a statement, not all of the subsequent 
works. Literature citation format: Authors (last name, 
followed by initials for first name and, if given, mid- 
dle name of first author; then list names of additional 
authors with initials before last names). Year. Title of 
article. Abbreviated title of the journal in which it was 
published. Always include either the range of page num- 
bers (for a journal article) or a total number of pages (for 
a book or other type of publication). List a sequence of 
citations in the general text chronologically, for example, 
“(Smith, 1932: Green. 1947; Smith and Jones, 1985).” 


Acknowledgments should be no more than 6 lines of 
text. Only those who have contributed in an outstanding 
way should be acknowledged by name. For recognition of 
other persons or groups, use a general term, such as crew, 
observers, or research coordinators, and do not include 
names with these terms. 


Digital object identifier (doi) code ensures that a publica- 
tion has a permanent location online. A doi link (which 
may include a doi code) should be included at the end of 
citations of published literature. Authors are responsi- 
ble for submitting accurate doi links. Faulty links will be 
deleted at the page-proof stage. 


Footnotes are used for all documents that have not been 
formally peer reviewed and for observations and personal 
communications, but these types of references should be 
cited sparingly in manuscripts submitted to the journal. 


All reference documents, administrative reports, inter- 
nal reports, progress reports, project reports, contract 
reports, personal observations, personal communications, 
unpublished data, manuscripts in review, and council meet- 
ing notes are footnoted in 10-point font and placed at the 
bottom of the page on which they are first cited. Footnote 
format is the same as that for formal literature citations. A 
link to the online source (e.g., [Available from http://www... , 
accessed July 2017.]), or the mailing address of the agency 
or department holding the document, should be provided 
so that readers may obtain a copy of the document. 


Tables are often overused in scientific papers; it is seldom 
necessary to present all the data associated with a study. 
Tables should not be excessive in size and must be cited 
in numerical order in the text. Headings should be short 
but ample enough to allow the table to be intelligible on 
its own. 

All abbreviations and unusual symbols must be 
explained in the table legend. Other incidental com- 
ments may be footnoted with numeral footnote markers. 
Use asterisks only to indicate significance in statistical 
data. Do not put a table legend on a page separate from 
the table; place the legend above the table. Do not submit 
tables in photo mode. 


e Note probability with a capital, italic P. 


e Provide a zero before all decimal points for values less 
than one (e.g., 0.07). 


e Round all values to 2 decimal points. 


e Use a comma in numbers of 5 digits or more (e.g., 
13,000 but 3000). 


Figures must be cited in numerical order in the text. 
Graphics should aid in the comprehension of the text, but 
they should be limited to presenting patterns rather than 
raw data. The number of figures should not exceed 1 figure 
for every 4 pages of text. 

Figure legends should explain all symbols and abbrevi- 
ations seen in the figure and should be double spaced on a 
separate page at the end of the manuscript. 

Line art and halftone figures should be saved at res- 
olutions >600 dots per inch (dpi) and >300 dpi, respec- 
tively. Color is allowed in figures to show morphological 
differences among species (for species identification), to 
show stain reactions, to show gradations (such as those of 
temperature and salinity within maps), and to distinguish 
between numerous lines and symbols in graphs. Figures 
approved for color should be saved in CMYK format. 

All figures must be submitted as PDF, TIFF, or EPS 
files. 


e¢ Capitalize only the first letter of the first word and 
proper nouns in all labels within figures. 


¢ Do not use overly large font sizes for labels in maps and 
for axis labels in graphs. 
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e Use the same point size for all labels, except for panel 
labels (e.g., A), which should be slightly larger than 
other labels (e.g., 11 point versus 8 point). 


e Use a sans serif font for all labels. Panel labels (e.g., A), 
however, should be in Times New Roman font. 


e Do not use bold fonts or bold lines in figures. Do not 
use italic fonts. Exceptions include use of italic fonts 
for labels of bodies of water in maps and a bold font for 
panel labels (e.g., A). 


e Do not place outline rules around graphs. 


¢ Do not include vertical and horizontal lines in the back- 
ground of graphs. Ticks for values on the x-axis and 
y-axis will suffice. 


e Place a north arrow and label degrees latitude and lon- 
gitude (e.g., 170°E) in all maps. If scale of map requires 
more than degrees, use degrees minutes, not decimal 
degrees. 


e Place panel labels (e.g., A, B, C) within the upper-left 
area of each graph or photo in a multi-panel figure, 
from left to right, then top to bottom. If the letter is 
not visible against a dark background, put a white box 
behind it. Do not use white labels. 


e Avoid placing labels vertically or diagonally. Y-axis 
labels can be vertical. Words in horizontal labels can be 
stacked vertically to fit. 


e Use symbols, shadings, or patterns (not clip art) in 
maps and graphs. 


e For scale bars in maps, use kilometers. Use the label 
km or kilometers (lowercased). 


Supplementary materials that are considered essential, 
but are too large or impractical for inclusion in a paper 
(e.g., metadata, figures, tables, videos, and websites), may 
be provided at the end of an article. These materials are 
subject to the editorial standards of the journal. A URL to 
the supplementary material and a brief explanation for 
including such material should be sent at the time of ini- 
tial submission of the paper to the journal. 


e Metadata, figures, and tables should be submitted in 
standard digital format (MS Word or PDF file) and 
should be cited in the general text, for example, as 
“".. was determined (Suppl. Table 3, Suppl. Fig. 1).” 


e Websites should be cited with a URL in the general 
text. 


e Videos must not be larger than 30 MB to allow a swift 
technical response for viewing the video. Authors 
should consider whether a short video uniquely cap- 
tures what text alone cannot capture for the under- 
standing of a process or behavior under examination 
in the article. Supply an online link to the location of 
the video. 


Copyright law does not apply to Fishery Bulletin, which 
falls within the public domain. However, if an author 
reproduces any part of an article from Fishery Bulletin, 
reference to source is considered correct form (e.g., Source: 
Fish. Bull. 97:105). 


Failure to follow these guidelines 
and failure to correspond with editors 
in a timely manner will delay 
publication of a manuscript. 


Submission of manuscript 


Submit a manuscript online at the ScholarOne Manu- 
scripts website for Fishery Bulletin. Commerce Depart- 
ment authors must provide proof of internal clearance 
of their manuscripts with either a completed and signed 
NOAA Form 25-700 or a copy of the clearance email from 
the Research Publication Tracking System. For further 
details on electronic submission, please contact the asso- 
ciate editor, Cara Mayo, at 


cara.mayo@noaa.gov. 


When requested, the text and tables should be submitted 
in MS Word format. Each figure should be sent as a sep- 
arate PDF, TIFF, or EPS file. Send a copy of a figure in 
the original software if conversion to any of these formats 
yields a degraded version of the figure. 


Questions? If you have questions regarding these guide- 
lines, please contact the managing editor, Kathryn 
Dennis, at 


kathryn.dennis@noaa. gov. 


Questions regarding manuscripts under review should be 
addressed to Associate Editor Cara Mayo. 
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| _] VISA or MasterCard Account your orders 
LOUIE eee eee Benepe 
[ | [| |] (Credit card expiration date) 


(Authorizing Signature) 


Mail to: Superintendent of Documents Thank you for 
P.O. Box 979050, St. Louis, Missouri 63197 oumonder. 
Tel. no. (U.S. callers): 866-512-1800 
Tel. no. (International callers): 1-202-512-1800 


Also available online at 
https://bookstore.gpo.gov/products/fishery-bulletin 
email: ContactCenter@gpo.gov 
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